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ABSTRACT 

We present a new method to calculate formation of cosmological structure in the New- 
tonian limit. The method is based on Lagrangian perturbation theory plus two key 
theoretical extensions. One advance involves fixing a previously ignored gauge-like de- 
gree of freedom present in the formal Lagrangian perturbation theory. The traditional 
derivation of the perturbation expansion introduces this unwanted freedom which it is 
crucial to eliminate. In effect, we transform the usual results of a perturbation theory 
calculation by a frame shift to give answers sought by a particular observer. A second 
extension is based on our previous work where we showed that Lagrangian pertur- 
bation theory expansions converge only over a limited time interval. The end of this 
interval generally precedes the onset of orbit crossing. We had introduced the idea of a 
multi-step method to extend the solution as far forward in time as possible. Here, we 
implement both the frame shift and the multi-step method to produce an algorithm 
capable of solving for the cosmological evolution of cold matter. 

Extensive 'proof of principle' tests validate the ideas and implementation of the 
method. These include direct comparison to known solutions, evaluation of conserved 
quantities and formal convergence studies (decrease of solution differences as the con- 
trol parameters are refined). The algorithm behaves satisfactorily in all these trials. 
The rate of convergence is exponential in the grid size, exponential in the Lagrangian 
order and polynomial in the step size. 

There are three main advantages afforded by the new technique. First and fore- 
most, it employs a smooth representation of all fields and the results are not limited 
by particle induced shot-noise errors. Second, the numerical error for any problem 
can be controlled by changing Lagrangian order and/or number of steps. In principle, 
arbitrarily small errors can be achieved prior to orbit crossing. Third, the initial data 
is completely generic, including cases where the initial velocity field has a rotational 
component. Together, these properties make the new technique well-suited to handle 
problems on quasi-linear scales where analytic methods and/or numerical simulations 
fail to provide suitably accurate answers. 

Key words: cosmology: theory - large-scale structure of Universe. 



1 INTRODUCTION 

An important means of determining and/or constraining cosmological parameters is to compare the growth of structure as 
predicted by theoretical models with that observed in the actual universe. The complexity and sophistication of theoretical 
approaches span an enormous range of possibilities. Linear perturbation theory provides the simplest analytic description of 
growth, generally suitable when the perturbation expansion parameter is small. Quasi-linear theory extends the accuracy and 
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applicability of the description at the cost of increasingly complex formulations. Finally, an N-body simulation is capable of 
handling fully non-linear situations and delivering, in principle, exact answers albeit at substantial computational cost and 
without the benefit of analytic insight. The choice of a particular method is best determined by the sort of problem one 
intends to solve. For example, if one is interested in calculating how an initial Gaussian perturbation spectrum generates 
non-Gaussianity as structure forms then one seeks an accurate method for treating perturbations once mode-mode coupling 
becomes important. As in Goldilocks's fable the extreme choices may prove unpalatable: on the one hand, the coupling 
of interest is absent in a purely linear analysis; on the other, particle-based methodologies introduce spurious shot noise 
in representations of the underlying smooth density distribution so that high accuracy may demand infeasible numbers of 
particles. 

This paper's focus is the development of a middle approach, a method for the numerical treatment of structure formation 
based on intrinsically smooth descriptions of the density and velocity fields and capable of tracing quasi-linear evolution to 
arbitrary accuracy. We anticipate it will find application in situations where small non- linear effects and/or highly accurate 
results prior to collapse are of primary interest. The physical context is Newtonian cosmology i.e. sub-horizon scales, non- 
relativistic velocities and weak gravitational fields. The methods are applicable until the formation of the first caustic. Future 
studies will extend the methodology to reach beyond particle crossing and to handle physical scales approaching that of the 
horizon. 

Lagrangian perturbation theory (LPT) , the heart of the approach employed here, has been well-studied. Its use in cosmol- 
ogy was initiated by Zel'dovich who gave an approximate analytic solution for special initial conditions with peculiar velocity 
proportional to peculiar acceleration (Zel'dovich 1970). This is the famous "Zel'dovich approximation" for the evolution of 
growing modes where "approximate" means accurate to first order in particle displacements. Buchert's studies ( Buchert|1992 



henceforth B92) enlarged the class of initial conditions that could be solved to first-order. Generalizations of the equations of 
motion, the initial conditions, the perturbation order and the method of solution have appeared regularly ( jMoutarde et al 
1991] IBouchet et al.||1992| I Buchert fc Ehlers|[l993l henceforth BE93; |Buchert|[l994l henceforth B94; |Catelan||l995l |Munshi 



Sahni, fc S tarobinsky 1994, Bouc het et aTp 995 Ehl ers fe Buchert|1997| henceforth EB97; [Rarnpf fc Buchert|2012| [Rampf fc 
Rigopoulos 2012] |Rampf fc W ong 2012] |Tatekawa|2012|. Physical extensions include general relativistic equations of motion 



( Kasai| 19*95 1 and the treatment of relativistic fluids ( Matarrese fc Terranova|[l996 Matarrese, Pantano, fc Saez|1993 1994 1 

LPT has proven to be a practical tool in many different applications including modelling the non-linear halo mass functions 
( Monacoj 1997 Scoccimarro fc Sheth|2002 \ , reconstructing baryon acoustic oscillations ( jEisenstein et al.|2007 \ , modelling the 
non-linear density velocity relation ( Kitaun^e^jflJ|2012 \ and setting up the initial conditions for fully numerical simulations 
( Scoccimarro||1998 Crocce, Pueblas, fc Scoccimarro|2006 1. 

An assumption common to all these treatments is that the matter velocity at a point in space is single valued (the 
distribution is "cold"). LPT breaks down once particle orbits begin to cross, caustics appear, densities diverge and physical 
singularities first form. Nothing in this paper ameliorates the onset of these effects. There have been interesting suggestions 



on how one might modify LPT to handle this fundamental limitation (e.g. Adler & Buchert 1999 1 but we do not pursue 
the issue here. For our purposes the immediate concern is a new limitation of LPT that has recently come to light. The 



convergence of the formal series expansion which is the basis of LPT turns out to be limited in time (Nadkarni- Ghosh & 



|Chernoff||2011| hereafter NC11). No one would be surprised that the onset of particle crossings inhibits the utility of the 
formal perturbation expansion for subsequent times. What was unexpected and striking to us was the existence of limits in 
a variety of other situations, even cosmologies lacking future particle crossings and devoid of future physical singularities. 



To the best of our knowledge Sahni & Shandarin ( 1996 1 were the first to report the failure of LPT to describe one of the 



simplest cosmological problems, the evolution of spherical homogeneous voids. We confirmed their result, investigated the 
issue for all top-hat cosmologies and traced the problem to the occurrence of poles in a suitable set of analytically continued 
equations of motion. We complemented this mathematical explanation with direct numerical evaluation of the LPT series in 
all qualitatively distinct regimes for top-hat like cosmologies. The upshot is the existence of a non-trivial limit to the future 
times for which the LPT description converges. We dubbed the interval when convergence of the perturbation expansion 
occurs as the "time of validity." Even if one could work to infinite order in the perturbation amplitude the LPT series would 
fail to give the correct answer at times outside the interval. The end of the time of validity precedes any limitation due to 
particle crossing. Figure [I] shows the first 15 orders of the LPT description for a simple underdense top-hat. The time of 
validity extends only up until a ~ 0.2. 

Like us many readers may find it odd that this behaviour had not been previously reported in a model as well-studied 
as the top-hat but we can point to two reasons. First, convergence of LPT had never been addressed in a systematic fashion. 
Second, for collapsing Zel'dovich initial conditions (initial velocity and acceleration proportional) the situation is degenerate 
in the following sense: the convergence limitations become identical to the moment of caustic formation. Generically, however, 
the time of validity for top-hat evolution does not extend up to the first particle crossings and is the more stringent restriction. 
To calculate the "answer" beyond the time of validity we proposed to evaluate the solution at an intermediate time when the 
series is still valid and use the results as initial conditions for a new expansion about the intermediate point. We worked out 
the details for this solution technique which we dubbed "LPT re-expansion" . It is roughly analogous to analytic continuation 
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of a power series in the complex plane. Re-expansion will not circumvent the ultimate limitation set by particle crossing but 
it will allow the solution to be extended to this maximum possible future time after which the flow is no longer cold. We 
provided detailed numerical examples of how LPT re-expansion handles the previous difficulties encountered in the top-hat 
problem. 

Although our analysis was specific to one of the simplest of all cosmological problems we believe these ideas are valid 
in more general circumstances. Consequently, one should shift from thinking of LPT as a one step, stand-alone method of 
calculation to viewing it as the fundamental operation in a numerical, multi-step method of solution. Our development of 
LPT re-expansion resembles a traditional finite-difference method in that the system is updated on a step-by-step basis. Just 
as is true for particle-based calculations, stability limits how big a step can be taken and accuracy varies with step size and 
expansion order. 

This paper works out LPT with re-expansion for general inhomogeneous initial conditions in a flat Friedmann-Robertson- 
Walker (FRW) background cosmology. We present fundamental mathematical results, practical computational methods and 
stringent numerical tests. This paper serves as a 'proof of principle' demonstration of the method. We will provide practical 
cosmological applications in subsequent papers. 

In S|2] we describe the Lagrangian framework and the relationship of the full solution of the geodesic equation of motion 
for particles to that provided by the LPT expansion. The formal expression for the latter is derived by first taking spatial 
derivatives of the geodesic equation (B92; BE93; B94; EB97). As we show in some detail, it is necessary to supplement the 
order-by-order solution of LPT to recover the complete solution of the geodesic equation. This issue that has never been 
addressed for single-step applications of LPT. For some initial conditions the omission turns out to be inconsequential but 
for most others it breaks momentum conservation. Moreover, when multiple steps are taken it becomes crucial to employ the 
full solution at each stage of the calculation. We will show that the numerical method converges if and only if that is done. 

This paper notably differs from most previous approaches by considering initial conditions of a completely general 
nature. LPT is often used to propagate the small fluctuations present at recombination to a modest redshift for the purpose of 
initialising an N-body treatment. In studies of cosmological structure formation it is common to begin from Zel'dovich initial 
conditions since vorticity modes decay away in an expanding universe. The damping is normally sufficient to obviate the 
influence of the initial vorticity so that an LPT treatment of the growing modes alone is adequate for setting up an N-body 
treatment of structure formation in pure ACDM. Vorticity is intrinsically generated once non-linearities form and particle 
crossing begins; these complications are accurately handled by traditional N-body methods. In this context, there is little 
need for treating arbitrary initial data. 

Of course having a capability to study both longitudinal and transverse motions opens up interesting possibilities involving 
additional participants in structure formation such as magnetic fields and cosmic strings since these may act to introduce 
vorticity into the cosmic flows at late times. But the main motivation to work in full generality is to facilitate the use of 
LPT as part of a future calculational approach meant to treat non-linear structure formation. In any realistic multi-step 
numerical approach (analogous to that developed herein) vorticity will inevitably appear just as it does in traditional N-body 
treatments. After it arises each step of a multi-step calculation must be able to handle essentially arbitrary density and velocity 
configurations as initial data. The ability to treat these initial conditions with LPT is a pre-requisite for utilizing LPT itself 
in the future approach. That is the main reason we work in full generality. 

Sj2]includes an outline of the design and implementation of the re-expansion algorithm, ^summarizes tests for both special 
and general initial conditions. We test the algorithm's ability to reproduce known solutions. We analyse the dependence of 
errors on Lagrangian order, step size and size of the numerical grid. Taken together these tests validate the concept of 
Lagrangian re-expansion, its analytic form and our computational implementation. All tests are for small grid sizes but we 
expect the methods demonstrated to apply without essential modification to larger scale calculations, 
presents conclusions and possible future developments. 



2 GRAVITATIONAL FIELD EQUATIONS IN THE LAGRANGIAN FRAMEWORK 
2.1 General Setup 

Assume a flat FRW background cosmology with scale factor a(t) that contains pressureless dark matter with average density 
pm{t) and time-dependent dark energy with average density px(t) having fixed equation of state parameter w. We do not 
distinguish between cold baryons and cold collisionless matter and we work before the formation of caustics. The homogeneous 
background evolution is determined by the initial scale factor ao, Hubble constant Ho and matter and dark energy densities, 
p m .o and px,o respectively. The subscript "0" indicates time to which may differ from the current epoch. 

We intend that the initial conditions allow for arbitrary density and velocity perturbations of cold dark matter subject to 
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a(t) 



Figure 1. Convergence of LPT: a(t) and b(t) are the scale factors of the background and the underdense spherical perturbation, 
respectively. The black dotted line is the exact open model with Sp/p = —9 X 10 -3 at a = 10 -3 and fractional Hubble velocity 3 X 10 — 3 . 
The blue lines show successive higher order approximations from 1 to 15 (as labeled). The LPT series converges only a<0.2 even though 
there is no future singularity in this cosmology. LPT re-expansion provides a method to extend the solution to arbitrarily large a in this 
example. 



minimal restrictions. We constrain the dark energy component to be spatially uniform at all times]]]. We work in an expanding, 
periodic volume, a 3-torus. For physical, cosmological applications this volume is intended to be (1) a fair representation of 
the actual universe and (2) one that can be accurately treated in the Newtonian limit. Our perturbative approach does not 
allow any back-reaction of the small-scale physics on the large-scale expansion of the homogeneous and isotropic background. 
Consistency requires that the initial mean density of the inhomogeneous matter distribution agree exactly with that of the 
background cosmology. Mathematically, the evolution of the inhomogeneous Newtonian problem is completely determined by 
the otherwise arbitrary, initial density and velocity perturbations. We must check, in principle, that the Newtonian limit is 
valid at all times. 



2.2 Lagrangian Picture 

Let us start in an inertial frame. In the Lagrangian picture the particle position is a function of time and a set of Lagrangian 
labels Y. The labels are assigned at the initial time and, for our purposes, are always taken to be initial comoving Eulerian 
positions. We write x(t; Y) for the instantaneous comoving coordinate of a fluid particle. By contrast, in the Eulerian picture 
x is an independent variable that identifies a fixed grid position. For notational simplicity we will often use the descriptive 
shorthand that x is either Lagrangian (equivalently, a particle) or Eulerian. 

In many applications one starts with a finite, discrete set of particles (e.g. an N-body code) but here we assume space is 
completely filled with particles. In effect, we promote the discrete set of Lagrangian labels to a continuous field throughout 
space so that at any time it is possible to relate an Eulerian grid coordinate to a Lagrangian label in a one-to-one fashion. 
Mathematically, the existence of a well-defined velocity field at every point in space establishes this result (EB97). The 
Lagrangian particle velocity is identical to the Eulerian velocity at the grid position that coincides with the particle position. 

Throughout this paper we will use the following convention for density: p(r, t) is the physical density expressed as a 
function of r, the physical position, and p(x, t) is the same numerical quantity expressed as a function of x, the comoving 
coordinate. Explicitly, we have 

p(x,t) = p(r,t)| r=ax (1) 
AM = p(r, t)d 3 r = p(x, t)d 3 r = p(x, t)a 3 d 3 x. (2) 

Note that p(x, t) is not the comoving density. 

Likewise, we employ a similar convention for density written in terms of the Lagrangian variable 

p(Y,t) = p(r,t)\ r=r (y t t) (3) 
AM = p{r,t)d 3 r = p(Y,t)d 3 r = p{Y ,t)J{Y ,t)d 3 Y (4) 



1 Although it has been suggested that quintessence models with w ^ — 1 should have spatial fluctuati ons (for example |Caldwell, Dave,| 
&: Steinhardt|1997 1 , these fluctuations are estimated to be small (for example Mota, Shaw, &: Silk|2008 Cooray, Holz, &: Caldwell|2010| . 
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J(Y,t) = Det(|^) (5) 

where r(Y, t) is the physical position of the particle labeled by Lagrangian coordinate Y and J is the determinant of the 
coordinate transformation. The distinguishing feature of the Lagrangian labels is that they track the mass so, in addition to 
the above general properties, we also have mass conservation in the form 

p(Y,t)J(Y,t)=p(Y,t')J(Y,t') (6) 
for any times t and t' . 

In comoving coordinates the Newtonian limit of the geodesic equation and of Poisson's equation for the gravitational 
potential ^>(x, t) are 



4- (« 2 ?) 

dt V dt J 



) = -V*V(M) (7) 

Viip(x,t) = 4nGa 2 Sp m (x.,t) (8) 



where a(t) is the scale factor and 5p m (x,t) is the physical matter density perturbation. There is also an equation for mass 
conservation. 

The a;-related terms are Eulerian except for those on the left-hand side of the first equation which refer to particles (d/dt 
and x). To make mathematical sense of these equations as a Lagrangian system we should, schematically, specify Lagrangian 
labels at the initial time and transform the Eulerian to Lagrangian forms while imposing the condition that the Eulerian 
coordinate match the particle position. We regard the system as posing an initial value problem for the particle positions and 
velocities. If the initial conditions are specified as Eulerian density and velocity perturbations then they must be transformed 
to particle positions and velocities. The density perturbation is a functional of the initial density and the initial and final 
particle positions. The potential can be solved at any time the density is known. 

In physical coordinates r = a(f)x we have 

av — fir = — aV r ?/>(r, t) (9) 
V?VM) = 47rG5p m (r,t), (10) 

where the dot denotes the total derivative with respect to time. As above, r(t) and r(t) (left-hand side of the first equation) 
refers to particles while the structure of the remaining terms is Eulerian. 

Buchert and Ehlers' derivation of the Lagrangian equations of motion (B92, BE93, B94, EB97) involves taking the 
Eulerian divergence and curl of the first equation, combining with the second and using the homogeneous equations for the 
scale factor a(t) to give 

V r -r = -4nG[p m {r,t) + p x (t){l + 3w)] (11) 
V r xi = 0, (12) 

where p m (r , i) and px (t) are the matter and dark energy densities respectively. These steps have the virtues of eliminating the 
potential and operator occurrences of x or r on the right-hand side. The essential complication, however, is that the derivative 
form of the geodesic equations is insensitive to a spatially homogeneous arbitrary, time-dependent vector shift Ar(i). 

The solution of the resultant Lagrangian system by Ehlers and Buchert (EB97) is achieved in a specific frame, one in 
which the initial, spatially averaged peculiar velocity vanishes. We must address the relationship between the solutions of 
the geodesic equation and of this specific solution of the derivative form of the geodesic equation. To handle this situation 
carefully, we introduce the observer frame (denoted by OBS) and the Ehlers-Buchert frame (denoted by EB). The solution to 
the geodesic equation lies in the former and the specific solution to the derivative form in the latter. The connection is called 
the frame shift. 

Schematically, the initial data is specified in the OBS frame, transformed to the EB frame, solved in the EB frame using 
EB's perturbative expansion and then transformed back to the OBS frame. 



2.3 Notation for Observer and EB Frames 

Begin by assuming that each frame possesses its own comoving coordinate system. Imagine a single particle variously de- 
scribed according to multiple coordinate systems. In each frame the Lagrangian label is assigned to be the comoving Eulerian 
coordinate at the initial time 

Y = x O bs,0 (13) 
X = x BS ,o. (14) 
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At future times the particles positions are written xobs(Y ,t) and x_bb(X, i), respectively. Physical coordinates, Jacobians 
relating Eulerian to Lagrangian coordinates and volume elements are defined in the usual manner in the OBS frame 

roB S (Y,t) = a{t)x Bs(Y,t) (15) 

J(Y,t) = Det(^) (16) 

d 3 x OBS = J -^fid A Y (17) 

and it follows J(Y,to) = ao 3 . Similar definitions are given in the EB frame. 

Without loss of generality choose the two frames to coincide at the initial time, xoss = "X-eb and Y = X, and write the 
relationship between the comoving coordinates as a pure function of time 

x.OBs(Y,t) = x BB (X,t) + Ax(f). (18) 
Here, Ax(t) is the comoving frame shift. Next, define Ar in terms of Ax 

Ar(t) = o(t)Ax(t). (19) 
We now have 

Ax(to) = Ar(t ( ,) = (20) 
J(Y,t) = J(X,t). (21) 

Initially, dPxoBS = d 3 Y — d 3 X — d 3 XEB, and at all times the Eulerian volume elements satisfy d?xoBS = d 3 XEB and the 
Lagrangian volume elements d 3 Y = d 3 X. We will choose Ax(io) in a moment and then later make use of the remaining 
freedom to choose Ax(t) so that the geodesic equation is solved. 

The physical mass density is pm{r,t). Since it's a scalar it remains unchanged by the EB to OBS transformation. Mass 
conservation and equations Q and (171 give 

p m (x.,t)a 3 d 3 x BS = Pm(Y,t )a 3 d 3 Y. (22) 
The peculiar velocity is defined 

voss(Y,*) = i Bs(Y,t)- a r OBS (Y,t) (23) 

a 

= a± Bs(Y,t) (24) 
and likewise for veb- The peculiar velocity shift is defined 

Av(t) = v Bs{Y,t) - v B fl(X,t) (25) 
= Ar(t) - -Ar(t) (26) 
= aAx(t). (27) 



2.4 Transformation from Observer to EB Frame 

The perturbation is characterized by two quantities specified in the observer frame; the initial fractional overdensity 
6oBs{Y,t )= pm{Y ' to) -1 (28) 

Pm.fi 

and the initial peculiar velocity 

voBs{Y,t ) = roBs(Y,t )-a Y. (29) 
The density perturbation is a scalar quantity so transferring to the EB frame is immediate 

5 ss (X,i ) = S O Bs(Y,t ) (30) 

and we will henceforth drop the "subscript" that distinguishes these two. 
The initial peculiar velocity is a vector quantity which transforms like 

v ES (X, t ) = v O Bs(Y,t ) - v Cj0 (31) 

where v Cj o is a constant we will now choose. Let (•••)« = (jd 3 s...)/Jd 3 s. We want to ensure that the average peculiar 
velocity in the EB frame vanishes at the beginning of the step, {veb)x = 0, as assumed in the formal perturbative scheme 
(see EB97). First define the average peculiar velocity in the OBS frame 
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v c (t) = (v Bs{x.OBS,t)) XOBS (32) 

and then note that at the beginning of a step the Eulerian and Lagrangian coordinate systems for both frames coincide, so 
that v Cj o = v c (to) will suffice. The choice implies that the initial time derivative of the frame shift is 

Ar(to) = aoAx(to) = v c ,o- (33) 
2.5 Equations in EB Frame 



Conservation of mass and equations 1 14 1, (281 and (30 1 imply 



(l + c?(X,fo))a 3 -J 

5(x '* ) = — (34) 

The density evolution of dark energy (denoted by a subscript X not to be confused with the Lagrangian coordinate) for a 
constant equation of state w is 

/an \ 3 ( 1 +™) 

px(t) = PXfi (^) ■ (35) 

These allow one to express the right-hand side of equation in terms of Lagrangian coordinates. We transform all the 



derivatives with respect to the Eulerian coordinates r = veb in equation (111 and (12 I to derivatives with respect to the 



Lagrangian coordinates (Appendix [A] furnishes additional details). The complete set of equations is 

L[r E B,r E B,r E B] = -3ff 2 ^m,oao 3 (l + S(X, to)) (36) 

2~(1 + 3w){l x ,o I — 1 L[t E b, veb, v E b\, 

T[r ES ,r BS ] = 0, (37) 



where 



InGp 



° m >o = orT^' . (38) 
3iio 

= 8^Gp^ 
3# 

and the action of L and T operators on general vectors A, B, C is defined as 

f,, „ „, dAi dBj dC k /(1 „s 

L[A,B,C] = eimq£ijk — ^L.— (40) 

T 9 [A,B] = (41) 

Here eijfc is the usual Levi-Civita symbol and Einstein's summation convention is used. The scalar operator L and vector 
operator T provide a compact representation. 

The specification of an, Ho, fim.o, ^x,o, w, the initial fractional overdensity and the initial peculiar velocity determines the 
solution. The fractional overdensity appears explicitly in the equations above; the peculiar velocity enters as initial conditions 
for r E B- 

2.6 EB Perturbation scheme 

In the formalism outlined by B92, BE93, B94 and EB97, the physical coordinate is split 

rs B (X,t) =o(t)X + p(X,t) (42) 

where the first term is the position in a homogeneous universe (zero density perturbations and zero peculiar velocities) and 
the second represents displacements that occur in the general case. The choice of the initial comoving coordinate as the 
Lagrangian label (equation |14[ | implies p(X, to) = 0. 
The displacement vector is expanded 

r EB (X,i) =a(t)X + ^p (n) (X,t)e n , (43) 

n=l 

where e is the formal, bookkeeping device that tracks the order of the displacement. Substitute the ansatz equation (|43[) into 



equations ( 36 1 and ( 37 \ and equate the terms of the same order of e to generate a hierarchy of equations to be solved. Since the 
Lagrangian labelling is independent of time the derivatives of interest commute: [d/dt, Vx] — 0. This makes it straightforward 
to find teb and veb- 
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This formal treatment at zeroth order of equation ( 36 1 reduces to 



" ff ° 2 ^^ + (l + 3»)«x,o(^) 3(1+ '" ) ) (44) 



a 2 



while equation (37 1 is identically zero. This is simply the equation governing the background scale factor. Given clq, Hq, Q. r , 



Q.x,o and w the background evolution is completely determined. 



At first order equations ( 36 1 and ( 37 1 reduce to 



^[Vx-p™] = -~H 2 n m , a 3 5(X,t ) (45) 

A T [V x xp (1) ] = (46) 
and at higher orders to 

Dt [V x ■ P (n) ] = S (n ' L) (47) 

D, T [V,xpW] = S<"' T \ (48) 
where the operators are 

A = I 2aa + -a H (1 + 3w)Ojc,o I — I d*a J ' ^ ' 

and S^" 6 ' and S (n ' T) are scalar and vector source terms comprised of combinations of displacements with order < n. An 
explicit form is given in the Appendix [B] Here L and T refer to longitudinal and transverse operators. 
Next, the displacement is split into longitudinal and transverse parts, order-by-order: 

In) (n.L) . (n,T) /ri\ 

p u = p v ' ' + ' (51) 

where V x x p( n ' L ' = and V x ■ p' n,T - ) = 0. On the 3-torus, a unique decomposition for order n requires (p' n ' (X.,t))x = 
(see Appendix C of EB97 for the mathematical proof). The choice of Lagrangian labels guarantees the condition is satisfied 
at t = to for all n. If the volume-averaged source terms vanish at all times then the decomposition is also possible at future 
times. Appendix [c] shows that (p^(X,t))x — at future times if both the overdensity and the peculiar velocity average to 
zero initially, which is the case (the latter by the choice of the EB frame). Using the n = 1 results one next shows that the 
average scalar and vector source terms for n = 2 displacements vanish so (p' 2 ' (X, t))x = and this argument is extended 
order-by-order. The result is that the requirements for decomposition are satisfied at all orders. 

Separating longitudinal and transverse displacements order-by-order uncouples the left-hand sides of equations (|47l) and 



( 48 1 and rationalizes the "L" and "T" labelling a posteriori. There are separate source term (S (n ' L) , S (n ' T) ) but note that 
each depends upon both types of lower order solutions (p( m,i ' and p( m ' T > for m < n). The entire solution at lower orders is 
needed to compute either the longitudinal or transverse source for a given order. 

The specification of the initial conditions for the hierarchy of equations determines the physical meaning of the formal 
perturbative expansion in powers of e. We make specific choices to ensure that the initial density and velocity perturbations 
are first order in e. The first choice is p( n ' L / T ) (X, to) = for all n ("L/T" means the equation applies to both longitudinal 
and transverse displacements) which is sufficient but not necessary to insure zero initial displacement. Homogeneous initial 
data makes it easy to use the formal structure of the hierarchy to check how powers of density and velocity enter. The density 
appears explicitly only at first order, i.e. in the equation for p^ 1 ' K The second choice is that the peculiar velocity is assigned 
as initial data only to the derivative of the first order term 

p (1 ' L/T) (X,t ) = v^/(X,t ) (52) 
where v^^(X, to) are the curl-free and divergence-less parts of the initial velocity. For all n > 1 

pK L/T) (X,t ) = (53) 

One might contemplate other choices (i.e., one could deviate from homogeneous for the displacements and/or distribute the 
peculiar velocity to different orders) but these would not have the simple physical interpretation that the natural choices 
provide and could vitiate the split into longitudinal and transverse components. 



In equations (451, (46 1, (47 1 and (481 the spatial and temporal operators commute and the spatial and temporal parts 
of the solution decouple at each order (see Appendix [5|. The solution is a sum of spatial parts times temporal parts. Each 
spatial part involves solving linear, elliptic partial differential equations of the form V ■ F = S L or V x F = S T with known 
sources S L , S T (Poisson-like problems). Each temporal part involves solving a second order ordinary differential equation in 
time (initial value problems). Both sets of equations include coefficients that depend upon the background cosmology. 
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As a practical matter all the Poisson-likc equations are solved using Fourier transforms on a TV x TV x TV grid with equally 
spaced grid points which represent the Lagrangian coordinates. Initial value equations are solved numerically using a standard 
differential equation solver. The individual displacement terms are summed to reconstruct r_es(X, t) using equation (43 1. 



2.7 Frame shifts 



The LPT expansion described in the previous section yields the position, velocity and density in the EB frame. To get the 
corresponding answers in the observer frame one must add the time dependent frame shifts (equations |18| and |25| l . 
Start with the Newtonian limit of the geodesic equation in the observer frame 



dt V dt ) 



j1p{X-OBS,t) 



(54) 



and take the average with respect to the Eulerian xoss (on the left we would formally replace dxoBs/dt with the local 
Eulerian velocity and d/dt with the full convective derivative). The right-hand side is periodic and averages to zero. Now 
formally we return to the Lagrangian description on the left: substitute for xoss from equation (18 1 and use d 3 xoBS = 
d 3 x E B = J(X.,t)a~ 3 d 3 X to give 



d_ 

dt 



,dAx(i) 
dt 



/ 2 dx E B \ 

V dt J 



a~ 3 J(X,t)d 3 X. 



and V — J d 3 x is the comoving volume of the box (physical volume a(t) 3 V). 

The geodesic equation is satisfied if the shift Ax satisfies this second order ODE with initial conditions Ax(io) 
Ax(to) = v Cl o/ao. The solution is 



Ax(t) = 

where 
g(t) = 




(t")dt" dt' + 



v c ,otso 



ar 



dt' 



dt 



(a x BS ) (a~ 3 J)d 3 X 



(op - Hp)(a- 3 J)d 3 X. 



Note that g(t) is completely determined from the EB solution. 
In the observer's frame the physical position and velocity are 

r Bs(Y,t) = aX + p(X,t) + aAx(i) 
r Bs(Y,t) = dX + p(X, t) + aAx(t) + aAx(t) 



and hence the peculiar velocity (defined by equation ( 23 ( ) is 



v Bs(Y,t) = p(X,t) 



-p(X,i) H h - 

a a a 



e(t')dt'. 



The density in the observer frame is identical to that in the EB frame 
5 Bs(Y,t) = 8 E B(X.,t) 

where the right-hand side is given by equation (|34[). 



(55) 
and 

(56) 

(57) 
(58) 



(59) 
(60) 



(61) 



(62) 



2.8 Frame shifts and momentum conservation 

Frame shifts ensure that the geodesic equation is properly satisfied. This is not just a mathematical requirement but a physical 
one. In this section we consider the implications for momentum conservation and drop the explicit OBS subscript for clarity. 

The total proper momentum in comoving coordinates written as a sum over mass elements AM, or an integral over 
Lagrangian labels (i.e. the density at the initial time) is 

i 

= f a^p m (Y,t„)agd 3 F. (64) 
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where x = a;(Y,t). Multiply by a and take the time derivative 
d(aP) 



dt 



Qp m (Y,t )a 3 d 3 Y, (65) 



2 dx 



Q = a' — . (66) 

Now replace the variable of integration Y with x — x(Y,t) at the time of interest. This recasts the integration from 
a Lagrangian to an Eulerian form. Using equation (22 I the densities are p m (Y, to)af ) d 3 Y = p m (x,t)a 3 d 3 x — (p m (t) + 



Sp m (x,t))a 3 d 3 x and in the last equality we have introduced the background and the fluctuating densities. The integral 
splits 

^p- = Pm {t)a 3 J Qd 3 x + a 3 J Q5p m (x,t)d 3 x. (67) 

When the geodesic equation equation |7| is satisfied Q is proportional to \7ip. Since ip is periodic the first volume integral 
with integrand proportional to Q vanishes. When the Poisson equation equation (JsJ> is satisfied one may rewrite using 
the formal Green's function for tp and invoke its symmetries to show that the second integral vanishes. Hence, d(aP)/dt = 
and P oc 1/a. 

If Poisson's equation is solved exactly in the EB frame the second integral will vanish. However, if frame shifts are ignored 
then the geodesic equation will not be satisfied and the first integral will generally fail to vanish. 

A qualitative understanding of why a non-zero v c must develop is revealed by looking at a problem with Zel'dovich 
initial conditions. We characterize this system as starting with first order terms in displacement, peculiar velocity and peculiar 
acceleration strictly proportional to one and other and negligible higher order terms (see Appendix|5]for a detailed discussion). 
Using v = ax, p m (t) = Pm,oCLo/a 3 and the definition of S, the total momentum is 

P(t) = Pm,oa 3 V (v c (t) + i J vS(x,t)d 3 x^J (68) 

where V is the comoving volume of the box. The initial state with velocity proportional to acceleration implies that the first 
term v c = because the system is periodic and the mean acceleration vanishes. The second term oc J Vip8d 3 x = on account 
of the symmetries of the Green's function for Poisson's equation. Hence, Zel'dovich initial conditions in a periodic system 
imply P(to) = 0. If/when the system evolves so that v is no longer proportional to acceleration then the second term will no 
longer vanish. The first term must be present and non-zero to insure P(t) = 0. In summary, the EB frame's non-zero velocity 
is intimately tied to the momentum conservation law which is maintained by proper inclusion of the frame shifts. 

Whenever a single step, first order scheme is used to initialize an N-body calculation starting with Zel'dovich initial 
conditions the frame shift vanishes, i.e. g(t) = 0, accurate at first order. To see this insert the definition of J from equation 



(A3 1, the proportionality of displacement, velocity and acceleration into equation (56 1. In general, a second order scheme that 
omits the frame shift completely makes a second order mistake, i.e. is accurate at first order but not second. 

This example also helps one understand some special situations when frame shifts are exactly zero. Assume the initial 
density and velocity field are purely one dimensional i.e., parallel planes of matter in 3D. As long as particle trajectories do 
not cross the acceleration on a particle is fixed by the initial density field. A single step first order LPT scheme yields an exact 
description and it is possible to write out explicit expressions for the frame shifts. For this one dimensional limit the mean 
velocity (in the OBS comoving coordinate system) is conserved i.e., {v) XOBS ~ v Ci oao/a. If it is initially zero, it is always 
zero. This does not imply that frame shifts vanish because there are two contributions, one coming from v c ,o and the other 



from internal dynamics in equation (561. In one dimension both terms will vanish if the initial velocity is proportional to the 
initial acceleration. In this special situation the observer and EB frames coincide at all times and the frame shift terms are 
exactly zero for both single and multi-step schemes. Proofs are provided in Appendix |E| 

Frame shifts are needed in other contexts and omitting them will generally produce inconsistencies that manifest as lack 
of momentum conservation. To the best of our knowledge these frame shifts have been ignored in most applications of LPT 
so far. 

Let's briefly consider the physical nature of v c . If the box contains a fair sample of the universe and if the comoving 
coordinate system coincides with the preferred FRW frame (one with an isotropic cosmic microwave background radiation) 
then one anticipates v c ,o = 0. On the other hand, if super-horizon motions are present (so the notion of a fair representation 
is not satisfied) and/or if the observer's frame is not equivalent to the preferred FRW frame then, generally speaking, v c .o is 
non-zero. 

Furthermore, in an inhomogeneous cosmology v c (t) is not proportional to the mean momentum in the box. Of course, 
the box's self-generated gravitational forces cannot impart net momentum to the box as a whole. If the initial net momentum 
is zero it will always be zero, however, this does not imply that v c will vanish on future Eulerian comoving grids if it vanishes 
initially. In particular, if one chooses v c (to) = then at later time ti one generally has v c (£i) 7^ 0. 

As a practical matter if one employs the EB method of solution one cannot circumvent the need to start each step by 
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finding the frame where v c vanishes. There is no single, adroit selection to be made at the calculation's start because the EB 
frame is not tied to a conserved quantity. 



2.9 LPT re-expansions 

The previous two sections describe the complete formalism for single step LPT; start with initial densities and velocities in 
the OBS frame, transform to the EB frame, solve for the displacements, move back to OBS frame and compute the densities 
and velocities at the end of the step. Two questions remain. How big a step can be made and what to do if that step doesn't 
encompass the total evolutionary time of interest? 

In NC11 we answered those questions in the context of spherically symmetric perturbations in an Q = 1 cosmology. We 
demonstrated that all convergence problems associated with the LPT series could be overcome by re-expanding the series 
in overlapping time domains with each domain subject to a time of validity criteria. The time of validity was determined 
rigorously as functions of the initial density and velocity perturbations. 

Here, we assume that the time of validity for the inhomogeneous evolution can be estimated by treating each point in 
the box as if it were an isolated top-hat. That is, we use the local density perturbation and local velocity perturbation to 
calculate the time interval that would be allowed for a top-hat. We find the maximum time that lies within all the individual 
intervals throughout the volume and set the time step accordingly. 

The fractional overdensity S and the fractional Hubble parameter 5 V are the specific parameters used to characterize the 
time of validity for the spherical perturbation. For generic inhomogeneous initial conditions we write the natural generalizations 

8 = *(Y,to) (69) 

5v = -^-V r -r(Y,t )-l (70) 

= ~W-v(Y,t ). (71) 

The values for 5 and S v are identical in OBS and EB frame. The time of validity T(S,5 V ) was determined in NC11; we use 
those results to find the minimum of T over the Lagrangian grid ^ 

A positive S at a point implies an overdense region and a positive 8 V implies an expanding region. From our experience 
with top-hat evolution we know that if the time of validity is set by an expanding region then the LPT re-expansion scheme 
has the ability to extend the solution arbitrarily far into the future. If it is set by a collapsing region then the LPT re-expansion 
can extend the solution to caustic formation. We remind the reader that the re-expansion scheme does not include any physics 
for treating hot (multi-stream) fluid. 

In LPT re-expansion the quantities calculated at the end of one step form initial conditions for the next step. Since each 
step begins with the Lagrangian coordinates equal to the comoving coordinates there is one additional computation that is 
necessary. A new Lagrangian coordinate Yi must be defined at the start of the next step; this is related to the Lagrangian 
coordinate of the previous step Yo by 

v r OBS (Yo,ti) 

Yl = ' (72) 



where ti is the starting time for the step and ross(Yo, £i) is given by equation ([59b. The final quantities from the previous 
step are known on a uniform grid in Yo; they are transformed so that the initial conditions are specified on a uniform grid 



in Yi . We solve equation ( 72 1 for Yo given uniform Yi and interpolate the final quantities at Yi . We refer to this step as 
regridding. 



3 NUMERICAL TESTS OF THE CODE 

This section checks and verifies the theoretical scheme outlined in ^2] by carrying out selective numerical calculations. The 
truncation error in the calculations depends upon the Lagrangian order (n), the number of time steps (Nt), and the size of 
the spatial grid (N s ). We will also refer to N t as the re-expansion frequency. 
The basic elements of the calculation are: 

(i) Use FFT (Fast Fourier Transform) methods to solve all spatial equations, 

(ii) Use Runge-kutta integration to solve the ordinary differential equations for the time-dependent coefficients (including 
the cosmology background), 

2 The spherical top-hat system has no transverse component so our treatment ignores that additional complication present in the 
inhomogeneous system. We see no evidence that this omission leads to any practical difficulty in stability or in convergence. 
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(iii) Use Simpson's rule for the quadratures needed to evaluate the frame shifts (eq. ( 56 1 ) , 

(iv) Use Newton's method for root finding and exact periodic interpolation during the regridding step (the most time- 
consuming task). 

Generally, we complete each individual task to machine precision. In the following section "interpolation error" refers to 
the net total of these four sources of numerical error. The total is typically dominated by the roundoff error in the solution 
of the spatial equations. These manipulate matrices with ~ N% elements. We observe residual errors at roughly the level 
^Machine where ^Machine is the machine precision. The interpolation error is small compared to truncation error in all 
results we will discuss. The truncation error is the main focus of investigation. 



3.1 Convergence and truncation error 

Consider evolution over a finite length of time prior to particle crossing in a periodic box of comoving size L. There are two 
qualitatively distinct LPT schemes that, in principle, are capable of converging to the exact answer: 

(i) Fix N t and increase n and N s without bounds. 

(ii) Fix n and increase Nt and N s without bounds. 

In both schemes it is necessary to refine the representation of all spatial functions (N s —¥ oo). In both it is necessary to respect 
the time of validity for evolution. The first scheme is satisfactory if Nt is sufficiently large; the second scheme which decreases 
the size of the time steps as Nt — > oo will eventually respect any non-zero time of validity condition. Hybrid forms in which 
Lagrangian order and number of steps increase together are also possible. 

Achieving a theoretical understanding of the scaling of the errors is important for two reasons. The numerical methods 
should demonstrate the expected scaling. Agreement between expected and actual scaling is valuable information when one is 
validating a computer code. In addition, a theoretical understanding helps make the most efficient choices of control parameters 
for carrying out practical calculations. 

In the past, LPT with Nt = 1 and fixed n has been applied to practical cosmological problems while the formal convergence 
issues have generally been ignored. Occasionally, a solution has been studied as a function of the Lagrangian order n for a 



single step as in the first scheme (Buchert et al.|1997l. There are numerous comparisons with N-body calculations (Buchert 



Melott, fc Weiss||1994| |Melott, Buchert, fc Weiss||1995| |Karakatsanis, Buchert, fc Melott||1997[ ) but these are not sufficiently 



accurate to address the issue of LPT convergence. We will provide a detailed look at convergence with respect to n and Nt. 

We treat N„ in a somewhat different manner. The spatial representation within the box should be refined as order and/or 
number of steps increase but we are unaware of any study of LPT that did so. Here, we will study one example with smooth 
initial conditions to establish that for sufficiently large grids the dominant error is set by Lagrangian order and is insensitive 
to grid size. For the remaining applications we argue (on a case by case basis) that grid-related errors are so small that 
convergence can be studied while holding N 3 large and fixed. 

This strategy is a reasonable but not rigorous approach to testing a method's convergence in a comoving, periodic volume. 
It is also important to note that making the comoving box larger may be crucially important for achieving accurate physical 
results. This is distinct from solving a given problem in a given box. For example, if the box holds a finite subset of a larger 
physical system (e.g. the universe) and one is interested in mean quantities of the larger system then one might need to consider 
ever bigger boxes. Or, if the physical system in the box obeys different boundary conditions (e.g. isolated not periodic) then 
one might need to increase the physical size. In these situations convergence to the physical answer of interest will require 
that the grid density N s /L and box size L increase without bounds. 



3.2 Initial conditions 

At the initial time to, the system is completely specified by the choice of cosmology, the initial density field <5oss(Y,to) an d 
the peculiar velocity voss(Y,to). Table [T] gives a list of various parameters for the test runs presented here. 

(i) All tests are done in an Q = 1 Einstein-deSitter cosmology over time intervals short enough that perturbations do not 
collapse to form multi-stream flows. 

(ii) The initial conditions satisfy v c ,o = 0. This choice of convenience allows us to characterize the peculiar velocity initial 
conditions completely and directly in terms of transverse and longitudinal content. It has no essential impact on the general 
conclusions derived from the testing. 

(iii) The tests cover both special and generic initial conditions. Special initial conditions possess non-trivial symmetries. For 
example, it is well-known that transverse modes decay with time in expanding cosmologies. Practical cosmological simulations 
typically begin by setting v|J s (X,to) = 0; such initial conditions are special by this criterion. We also consider generic initial 
conditions to exercise the full dynamics. 
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Table 1. Table of parameters for the numeric tests presented in *J3| Here, n is Lagrangian order, Nt is number of geometrically spaced 
steps, N s is spatial grid size, At is the total time, to is the initial time. The rms sub-scripted quantities are root mean square of 8, 
the initial density contrast, F$, the initial first order acceleration, and v( L / T ) , the initial velocity (longitudinal, transverse). Each test 
is discussed in the Section noted. ID means a one dimensional density and velocity profile. 'Top-hat' means a compensated top-hat 
smoothed by a fixed Gaussian and including a non-radial, smooth velocity perturbation. 'Exact' refers to the analytic solution given for 
special initial conditions (Model I in |Buchert et al.|1997[ l. 



Section Type 


n 


N t 


N s 


At 

to 


$rms 


-pS 

rms 


v L 

rms 


rms 




3.3 




ID 


1 


1-1 


16 


0.75 


0.007 


0.112 


0.007 







3.4 




Exact 


3 


1 


16 


.5 


1.22 


0.11 


0.11 







3.5 




Top-hat 


2 


1 


24-64 


0.75 


0.72 


0.05 


0.03 


0.042 




3.6.1 


Zel'dovich 


1-3 


1-5 


16 


0.75 


0.077 


0.687 


0.687 







3.6.2 


Generic 


1-3 


1-6 


16 


0.75 


0.077 


0.687 


0.270 


0.523 



3.3 Evolution of ID perturbations 



First order LPT is well-known to be an exact method for treating one dimensional problems prior to particle crossing ( Peacock 



1999). An "exact" solution is directly found by quadrature irrespective of the identification of the EB frame. Here we compare 
exact results of this sort to answers generated by our general 3D numerical implementation of LPT re-expansion. 

This comparison serves as a minimal check on the full code at first Lagrangian order because it involves all the basic 
elements of the most general calculation: the Fourier representation of spatial functions, solution of the ordinary differential 
equations for the time-dependent coefficients in the LPT expansion, evaluation of frame shifts and the interpolation from 
one Lagrangian grid to another. The solution of this simple problem exercises the entire re-expansion LPT methodology. No 
special modifications were made to the 3D code for carrying out this test. 

The initial density and velocity perturbations were taken to be one dimensional functions consisting of a single Fourier 
component <5(X,to) = 0.01 cos 2-kXi/L and v(X,to) = 0.01{cos 2-kXi/L, 0, 0}, where L is the comoving length of the box 
and X\ is the component along the i-axis. With this choice the initial velocity and acceleration are not proportional and we 
expect non-zero frame shift terms. Calculations were made at Lagrangian orders n = 1 — 3, for Nt = 1 — 4 steps covering a 



fixed total time At with time steps ti = to(l + At/toY^ Nt for ^ i ^ N t . Regridding was done by solving equation (72 I after 
each step. 

Exact and numerical results are compared in the observer's frame at the final time. As expected, the contributions to 
the solution at second and third order were zero to machine precision. The left panel of figure [5] shows the errors (maximum 
absolute velocity difference on the grid between exact and numerical results) as a function of the numbers of steps taken to 
reach the final time. The lower dotted line at roughly the round off level of computation (marked "with frame shifts" ) shows 
that essentially perfect agreement with the exact answer is achieved for both single and multiple steps. This result validates 
(at first order) the implementation of LPT re-expansion. Note that a sequence of exact first order calculations is exact. 

Frame shifts play an essential part in the numerical calculation. We repeated the calculation except that we imposed 
Ax = Ax = in equation (591 and equation (60 1. This omits the frame shifts. In left panel of figure [5] the upper dashed line 



shows the results obtained. There is an easily detected error with respect to the exact answer, an error present even for single 
step LPT. 

Another useful comparison may be made between numerical single and multi-step calculations even if one is ignorant 
of the "exact" solutions. As we have already shown above when frame shifts are included single and multi-step first order 
LPT calculations are equivalent. However, if one omits the frame shifts this equivalency is broken. The solid line shows the 
difference between the numerical results for runs of a single step and those of multiple steps when all frame shifts are omitted. 
The explanation is simple: none of the answers is exactly right and the discrepancies are tied to different numbers of missing 
frame shifts. 

The right hand panel of figure [2] shows a test of momentum conservation for calculations with and without frame shifts 
and for varying numbers of steps. With frame shifts the numerical calculation achieves round off levels of accuracy; without 
them the deficiency in the conservation law is easily detected and, moreover, does not improve as the interval in time is refined. 

Taken together these simple tests help validate the method and underline the crucial importance of including the frame 
shifts. 
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Figure 2. Left hand panel: the dotted line shows the maximum difference of numerical and quadrature-derived, exact answers when 
frame shifts are included in the numerical answer; the dashed line shows the same information when frame shifts are omitted from the 
numerical answer. Only when shifts are included do the numerical answers agree with the exact ones to machine precision. The solid 
line shows differences in single and multi-step numerical answers in the absence of frame shifts (none of these are "exact"). Right hand 
panel: test of momentum conservation for the numerical runs with (dotted) and without (dashed) frame shifts. 

3.4 Exact solutions 

An excellent check is Model I in |Buchert et al.j ( |1997[ ). It is an exact analytic result for one step in the EB frame at the first 
three orders for special initial conditions. Our numerical results agree with the analytic results to machine precision for all 
three orders. The special symmetry of the initial data implies that many of the spatial functions of the full general scheme 
vanish. The first, second and third orders comprise 2, 3 and 16 non-zero terms, respectively, and correspond to 2/3, 1/3 and 
1/4 of the total number of possible terms at each order. The agreement using the full set of possible terms confirms that the 
both the time and space pieces of non-zero terms are encoded correctly. The agreement provides a consistency check for the 
spatial terms that vanish but says nothing conclusive about the temporal parts of these particular terms. 

This check is essentially independent of the question of frame shifts. The analytic expressions for the frame shifts are 
exactly zero at all three orders. The numerically derived frame shift gave the same result to machine precision. It is not clear 
if the frame shifts continue to be zero at all higher orders, but for the purposes of this test we are concerned only with the 
first three orders. 

Ideally one would construct more general test cases to check all possible terms in the code. Unfortunately, the method 
of constructing 'local forms' (B94) is cumbersome and involves guessing the correct solution for the Lagrangian displacement 
at each order so that the divergence and curl match the source terms order-by-order. We have utilized the checks against all 
the published solutions of which we are aware. 

3.5 Top-hat with fixed smoothing 

The idealized configuration for this test is an overdense sphere surrounded by a compensating region and an infinite homoge- 
neous background. If the initial velocity perturbation is zero the solution is simple and analytic. There are two practical issues 
that potentially interfere with running tests starting from this sort of initial data. The computational volume is intrinsically 
periodic and the computational method is designed to handle continuous not discontinuous spatial distributions. 

A periodic configuration necessarily includes interactions between neighbouring boxes because exact spherical symmetry 
is not maintained in an approximate numerical solution. A suitably compact configuration limits these errors and the difference 
between the isolated and periodic versions can be made negligibly small. [^] In principle, the top-hat is an excellent problem for 
validating the LPT expansion since the (isolated) solutions for the displacement and density can be computed analytically order 
by order. Numerical results can be compared to analytic ones at a given order (NC11). The significant complication is that the 
exact top-hat profile is discontinuous along the boundaries between the overdense sphere, the vacuum compensating region and 
the outer region of average density. Such discontinuities induce Gibbs phenomenon in the spatial Fourier representation and 
mask differences between expected and observed solutions. In principle, one can make the grid sufficiently large to overcome 

3 The Fourier representation of the density guarantees that the mass within each box is exactly equal to the background value. The 
identical symmetries of the Cartesian grid within the box and of the arrangement of the surrounding periodic copies guarantee that the 
dipolc interactions between copies vanish. Therefore, the leading interaction that distinguishes the isolated from the periodic problem 
is a quadrupole term. It is ~ 10 _0162JVs /L 4 for the configuration studied. This is roughly equivalent to machine precision at N 3 = 100 
and is very small compared to the truncation errors studied here. 



© RAS, MNRAS 000, [T]-?? 



Non-linear evolution using LPT re-expansions 15 



CO - 6 - 

° 8 



Interpolation error 



15 20 25 30 35 40 

N s 

Figure 3. Convergence of the solution as the size of the spatial grid or the Lagrangian order is increased. Cauchy differences of density 
on the OBS Lagrangian grid are plotted. The dotted line traces the differences for calculations having increasing numbers of grid points 
(N s and N s + 8) at fixed Lagrangian orders (n = 1, 2 and 3). In fact, the calculated differences, labeled by n, all overlap. The straight 
dotted line indicates the exponential convergence expected for a spectral method. The solid, horizontal lines are Cauchy differences 
with respect to Lagrangian order (An refers to orders n and n + 1) at fixed N s . Both types of errors are generally present. For small 
grids the grid-related spatial errors dominate and for large grids the Lagrangian-related finite order errors dominate. The crossover from 
A^-dominated to n-dominated errors occurs at TV*. Additional calculations (not shown) demonstrate that N* is smaller for smoother 
initial data. All the errors shown are much larger than the interpolation error shown by the dashed-dotted line. 



this interference but we do not attempt to do so here because the convergence is very slow. That means we must rely on 
indirect checks that the higher order terms in the LPT expansion are correct. 

Here, we smooth the top-hat density field by a Gaussian with fixed width a and also impose a smooth velocity field. 
The resulting velocity profile has both longitudinal and transverse velocities. The details of the configuration can be found in 
Appendix [F] The exact solution is unknown so we discuss convergence using Cauchy differences for quantity / 

£ f (a;a')= y/((f' - /") 2 }„ (73) 

where a and a' label one or more of n, Nt or N s (with the other parameters held fixed) and q is the coordinate system used. 
We take / = 5 for density and v for velocity at the final time. For / = v all three velocity components are included. 

We calculate a single step of small size so that the final time is well before collapse. The assessment of convergence 
involves Cauchy differences formed with runs having different grid resolutions and Lagrangian orders, e.g. a — {N s ,n} and 
a' = {N 3 + 1, n + 1}. The spectral representation of spatial functions suggests that if grid errors dominate then the Cauchy 
differences will decrease exponentially with N s . Conversely, if errors in the series expansion in LPT dominate then the scaling 
will be exponential with n. This implies a sharp transition from errors dominated by grid resolution to errors dominated by 
Lagrangian order. 

Runs with 16 ^ N s ^ 48 and 1 ^ n ?C 3 were carried out for fixed smoothing. We evaluated the Cauchy differences as the 
grid was refined (£(N S ; N 3 + 8) at fixed n) and then as the Lagrangian order was increased (£ (n; n + 1) at fixed N 3 ) for q = Y 
(on the OBS Lagrangian grid). Figure[3]shows the results for the density. The Cauchy differences for varying grids are virtually 
identical for n = 1 — 3 (this is not obvious from the figure, because the lines for different n overlap). The Cauchy differences 
for varying Lagrangian do not depend upon N s (this is obvious from the fact the Cauchy differences are basically horizontal 
lines). Define the grid crossover N* (n) to be the value of N„ that satisfies £({n, N s }; {n, N 3 + 8}) = £ ({n, N„};{n + 1, N 3 }) 
for a given n. The behaviour of the formal measure of convergence is 

F(X vWnilWilrtl ( £({n,N s };{n,N s + S}) if N s < 7V a * 

The grid-related refinement dominates in the first case and order-related refinement in the second. The figure makes clear 
that the spatial differences decline exponentially with N s (as they should for a spectral method) and the analytic structure for 
the LPT expansion suggests that order errors decline exponentially with n (this will be tested more directly in a succeeding 
section). One surmises that for two sources of exponential errors, the transition N* should be a linear function of n. 

The specifics of the transition from grid to order-dominated convergence are problem dependent. We have investigated 
initial conditions which differ from those above only in the degree of smoothing used (choice of a). We find iVJ is larger when 
initial data is less smooth and vice- versa. 

The practical upshot is that if the solution is smooth and N s > Ns then the spatial errors are much smaller than changes 
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associated with incrementing the Lagrangian order or decreasing the step size over finite ranges of n and Nt, respectively. We 
attempt to conduct all the rest of our tests in the limit N 3 > N 3 when the grid errors are negligible. 

3.6 Convergence with Lagrangian order n and number of steps N t 

This section presents tests of convergence with Lagrangian order n and number of steps Nt for initial conditions arising from a 
random Gaussian field at fixed spatial grid size. Two types of initial data are considered: Zel'dovich initial conditions (velocity 
proportional to the acceleration field; irrotational flow) and generic initial conditions (initial velocity chosen arbitrarily from 
random Gaussian fields; transverse and longitudinal components present and of comparable size). We focus on the truncation 
error for Lagrangian order n and re-expansion frequency Nt- We choose grids, initial conditions and time intervals for evolution 
for which we can be reasonably sure that grid-related errors at fixed N„ are so small that they do not interfere with the Cauchy 
differences for varying n and Nt- 

We proceed as follows. The power in the initial data is limited to frequencies less than half the Nyquist frequency. In this 
discussion, we refer to wavenumbers as 'frequency' (k = 2n/X), not to be confused with the re-expansion frequency Nt. For 
box-length L and grid size N s , the magnitude of the Nyquist frequency is nN 3 /L. Since gravitational dynamics is intrinsically 
non-linear we know that the power in these initially zeroed modes will grow as collapse proceeds. Ultimately, any power that 
reaches or exceeds the Nyquist frequency will manifest as error since the Nyquist frequency is fixed because the grid is fixed. 
The size of the power in the vicinity of the Nyquist frequency is taken to be a surrogate for error incurred by using a finite 
grid. We monitor the power that builds up to make sure that the spatial errors remain negligible. 

The power in the Nyquist frequencies is 



V s w 

where (k) is the set of wave numbers having a Cartesian component equal to the Nyquist frequency (or, we can monitor the 
power within a given range of the Nyquist value). Consider, for example, the Cauchy difference with respect to N t on the 
OBS Lagrangian grid. This may be rewritten using Parseval's theorem ( Press et al.|2002 l 



£f(N t ;N t ')=(^^\f^'pq-f N *(X)A =(^Y / \f Nt '^)-f Nt (k)\ 2 ^j , (76) 

where / is the Fourier transform of / defined as f m = A^^X^i /,e 27ri(m - 1)(i - 1)/iVs . The maximum contribution to the 
total that comes from the Nyquist modes is 




^jV?£ l/Wt ' (k) " '""'Ml 2 * 2 « /V* > J/ V ' sk >i-' - (77) 

where the Nyquist components for Nt and Nt are estimated to be of the same order of magnitude. When £f(N t ;N t ) >> 
Pf,N yq we infer that the uncalculated finite grid errors are small compared to the calculated Cauchy differences. Likewise, 
when £f(n;n') » Pf,N yq we infer grid errors are small compared to the Cauchy differences for n. The fractional overdensity 
and peculiar velocity are not independent. We always use max/ = {a i „} Pf,N yq in the figures that follow. 

3.6.1 Zel'dovich initial conditions 

The initial density is chosen to be a specific realization of a random Gaussian field and the initial peculiar velocity is set 
to be proportional to the acceleration at each point (Appendix The initial conditions satisfy v c (to) = 0. We calculated 
evolution over the same time interval for Lagrangian orders n = 1 — 3 and frequency of re-expansion Nt = 1 — 5. 

Before presenting the convergence results let us first discuss an important feature of this flow. The velocity is irrotational 
at the initial time in the observer frame. Kelvin's theorem states that for an ideal fluid under the influence of a conservative 
force, the velocity circulation around a closed contour comoving with the fluid stays constant in time. As a consequence, the 



vorticity of the velocity field is conserved. In particular, if the fluid is irrotational initially it stays irrotational (Landau & 
Lifschitz|1998| EB97). 



Kelvin's theorem applies to the exact dynamics in the Eulerian coordinate system. Although the initial conditions are 
irrotational in both Eulerian and Lagrangian systems, the Lagrangian calculations do not preserve irrotational flows in the 
Lagrangian coordinate system because there is no Lagrangian version of Kelvin's theorem. If the Lagrangian calculation 
were exact then the flow would satisfy the exact conservation law for circulation after transformation back to the Eulerian 
coordinate system. 

Our calculation is a finite order Lagrangian calculation with finite step size. We expect the exact Eulerian result that the 



© RAS, MNRAS 000, [fj-?? 



Non-linear evolution using LPT re-expansions 17 



longitudinal 





' transverse 



q 1 ~ 10 X Nyquist error 



constant 3 , 



Interpolation error 



Log 10 W t 

Figure 4. Root mean square of v' L ' and v' T ) and v c at the final time in the observer frame for the Zel'dovich initial conditions as a 
function of Nt the number of steps taken. Results for Lagrangian orders n = 1 — 3 have different colours and plot markers. Note that 
transverse components converge to zero in agreement with Kelvin's theorem. 



flow should be irrotational (at any time) will emerge only as convergence is achieved. All non-zero transverse components in 
the OBS frame are errors and should necessarily diminish as the calculation is refined. 

Figure [4] summarizes the scale of different velocity components at the final time in the OBS frame as a function of Nt for 
calculations at Lagrangian orders n = 1 — 3. In all cases the longitudinal component dominates and the transverse term is 
small by comparison. The transverse term decreases as Nt increases and as n increases. Both these trends show that solution 
behaves qualitatively in accord with Kelvin's theorem. 

Three other important scales are indicated on figure [4] The interpolation error is very small, consistent with the effect 
of machine precision; it is inconsequential for all our considerations. The Nyquist error provides an estimate of the error 
associated with the spatial grid. It is considerably larger than interpolation error but less than the scale of the transverse 
velocities mentioned above. Finally, the mean peculiar velocity v c (labeled "constant" because it does not depend upon 



position) is shown. As discussed in [ 2.8 v c is generally non-zero. In this run, the Zel'dovich initial conditions dictate v Ci n = 
and the time interval of evolution is small. Hence, the final v c is small in the sense of being no bigger than the Nyquist error. 
We infer that frame shifts play no essential role in these particular results. 

Now consider the Cauchy differences for the solutions in figure [5] The upper panel presents the results for density, 
longitudinal and transverse components of velocity (left to right) in the OBS Eulerian frame as the number of steps Nt 
changes. Lines are labeled by Lagrangian order n. The Nyquist and interpolation errors lie below the Cauchy differences 
indicating that the truncation errors for order and frequency of re-expansion have been isolated from the effects of grid size 
and numerical precision. The slope of the lines slightly increases with higher n indicating faster convergence. 

For the top-hat it was found theoretically, and confirmed empirically, that for a small time step At /to << 1 the Cauchy 
differences decreased with Nt and n as 

£to P -hat(N t ; N t + 1) ~ |ffjv t +i,n — QN t ,n\ (78) 
Stop-hat (n; n + 1) ~ \gN t ,n\, (79) 



where 

•ix „ = K Nt>n co*0 su\" OA" ' ' ' 



" +I (f?r < 8 »» 



KN t ,n = 

The quantities A and 6 characterize the strengths of the initial density and velocity perturbations (8 = Acos# and S v = 
A sin (9). Note that <?jv t ,n scales exponentially with Lagrangian order n for any given Nt. This is expected because the n-th 
term in the LPT series is oc e™, where e is the magnitude of the initial perturbation. Conversely, the theoretical error scales 
as a power of the number of steps oc iV t _n as Nt increases. The latter also implies that the slope steepens with Lagrangian 
order n. 

The lower panel of figure [5] presents a quantitative comparison of the top-hat convergence rates to the more general 
problem that has been calculated. The figure shows the ratio of Z to that deduced previously for the top-hat. (For Zel'dovich 
initial conditions 8 — arctan(— 1/3) and A = VT05/3; we used the root mean square value for 5). The ratio is approximately 
constant implying that the convergence rates with respect to n and Nt for inhomogeneous initial conditions are close to those 
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Figure 5. Convergence with respect to frequency of re-expansion for Zel'dovich initial conditions. Top panel: Cauchy differences at the 
final time £ v l, £ v t (comparing N t , N t + 1 for fractional overdensity, longitudinal peculiar velocity and transverse peculiar velocity, 
respectively). Nyquist error and interpolation error levels lie below those of the Cauchy differences and do not interfere. For all variables, 
at any order, the error decreases as the Nt increases. The rate of convergence is better for larger Lagrangian order n. Bottom panel: the 
ratio of Cauchy differences in the inhomogeneous calculation to those of the top-hat. It is noteworthy that the ratio is roughly constant 
for larger Nt and n. This provides a confirmation that the theoretical top-hat convergence rate describes the inhomogeneous case. 



of the top-hat. The agreement is very good for the density and longitudinal velocity. The only anomaly appears to be related 
to the relative size of the first and second order transverse velocity differences (far right panel) . The second and third order 
terms scale as predicted theoretically. 

The explanation is enlightening. One might first imagine that the cause was related to the fact that the top-hat collapse 
does not include any components of transverse velocity. Then it would be unsurprising if the theoretical result failed to capture 
the apparent inversion in transverse components. But that is not what we find. The key point is already present in figure [4] 
which shows that the magnitude of the transverse velocity is very small after a single, first order step but much larger after 
a single, second order step or after multiple steps of first order. 

We find the explanation is related to the way in which the Lagrangian system "recovers" the result implied by Kelvin's 
theorem, i.e. that the transverse component in the Eulerian system should ultimately vanish. When the initial conditions 
satisfy velocity proportional to acceleration, the lowest order, non-vanishing contribution to the Lagrangian frame transverse 
velocity is third order (BE93; B94). However, the transformation from Lagrangian to Eulerian coordinates is essential, i.e. 
the transformation itself can contribute to the Eulerian transverse velocity. The Eulerian transverse velocity (defined as 
w_E(r,t) = V r x voss) is equal to the Lagrangian transverse velocity (defined as wi(Y,f) = Vy x voss) plus 'extra' terms 
that are quadratic and cubic combinations of the displacement vector and its time derivative (see Appendix [G| 

w E (Y,t) =w L (Y,t) + Terms[p,p] + Terms[p, p, p]. (82) 

When the calculation is done with a first order, single step LPT scheme for Zel'dovich initial conditions, these 'extra' terms 
vanish exactly and we get both Wi(Y,i) = and Wb(Y,4) = at all times. When the calculation is done with a second 
order, single step LPT scheme for Zel'dovich initial conditions, wi(Y,f) = at all times, but the 'extra' non-zero terms are 
third and higher order in the expansion parameter ('extra' terms that are second order vanish). 

w E (Y,t) =w L (Y,t) + 0(e 3 ) (83) 

w_e(Y,i 1 ) 7^ except for the initial time. So, we have the unexpected consequence that a single step of a first order LPT 
calculation will yield zero Eulerian transverse velocity in the OBS frame but that a single step of a higher order calculation 
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Figure 6. Root mean square of various components of the final velocity in the observer frame (left) and the transverse velocity as a 
function of the scale factor (right). The longitudinal and transverse components are roughly the same magnitude initially and, on account 
of the small time advanced, finally. However, on an expanded scale magnitude of the transverse velocity decreases oc 1/a. Dots are data, 
line is the fit. The theoretically expected slope is —1 and best fit line has slope —0.99. All orders have similar behaviour; first order data 
is shown here. 



will generally yield a non-zero Eulerian transverse velocity. Both the first order and the infinite order calculations will satisfy 
Kelvin's theorem. Of course, the first order solution is not exact and we can expect convergence only as the order increases. 
In general, with a n-th order calculation, one expects w_e(r, t) — 0(e n+1 ) which relies on complicated cancellations between 
all terms of order less than n between the three pieces on the r.h.s. of equation (82 1. 

Figure [4] shows that multiple steps of a first order scheme builds up a non-zero transverse velocity. This might appear to 
be at odds with the discussion above in which a single, first order step exactly satisfies Kelvin's theorem. The reason is that 
even one step will generate a state where the initial proportionality between velocity and acceleration is broken. In subsequent 
steps even a first order LPT scheme can then give non-vanishing contribution to transverse velocity in the Eulerian frame 
(due to the second and higher order terms in the transformation between Eulerian and Lagrangian coordinates). Multi-step, 
first order calculations may generate transverse velocity. Convergence and consistency with Kelvin's theorem is expected as 
the number of steps increase. 

Ultimately, only the asymptotics matter: we have shown that refinements in order and/or step size both lead to smaller 
transverse motions in the Eulerian frame. Higher order Lagrangian schemes converge more quickly in the asymptotic limit. 
For example, in the top right hand panel the line for n — 3 lies below that for n = 2 and its downward slope (as a function 
of N t ) is greater. 

We conclude that the LPT re-expansion converges for Zel'dovich initial conditions and, moreover, that it does so at a 
rate which is essentially identical to that theoretically derived in the case of the top-hat. 



3.6.2 Generic initial conditions 

In this test the initial density and peculiar velocity are independent realizations of a Gaussian random field. The latter has 
longitudinal and transverse components of comparable magnitude. We follow the evolution of the initial conditions for varying 
Lagrangian order (n = 1 — 3) and re-expansion frequency (Nt = 1 — 6). Figure [fi] summarizes (left panel) the root mean square 
velocity at the final time and (right panel) the evolution of the transverse velocity. For the small time step considered here the 
longitudinal and transverse parts are still comparable at the final time but (right panel) the transverse mode decays oc 1/a. 
Here, v c is greater than in the Zel'dovich case and points to the importance of frame shifts. 

Figure [7] shows the Cauchy differences for changes in Nt . The differences decrease with increase in Nt and Lagrangian 
order n. The transverse term here behaves in a qualitatively identical manner to the longitudinal case (there is no low-order 
anomaly like that seen in the Zel'dovich case). The scaling of these errors closely follows the trend established for the top-hat 
(the graphs are not repeated here). 
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Figure 7. General initial conditions: convergence with respect to Lagrangian order and frequency of re-expansions. The coding of the 
graph is the same as figure [5] Cauchy differences for Aft decrease as Nt increases; the rate of convergence improves for higher Lagrangian 
order n. The scaling of these errors with respect to a top-hat is similar to the lower panel of figure [5] and is not shown separately here. 
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Figure 8. Frame shifts and momentum conservation (general case). Top panel: Cauchy differences for density with (left) and without 
(right) frame shifts. Bottom panel: momentum conservation errors with (left) and without (right) frame shifts. Frame shifts are crucial 
to ensure convergence with Lagrangian order and ignoring them leads to a breakdown of momentum conservation. 



Our main conclusion drawn from the results for the generic test case is that the full gravitational dynamics has been 
successfully implemented in the LPT re-expansion method. 



3.6.3 Frame shifts again 

We had previously illustrated frame shifts in the simplest possible case, one-dimensional problems. Here we revisit the issue 
in the context of the general case just discussed. We will show that the numerical convergence as measured by the Cauchy 
differences is destroyed and momentum conservation broken when frame shifts are omitted. 

Figures [8] compares convergence of the density field in calculations with (left) and without (right) frame shifts. The left 
panel shows the same sort of scaling with n and Nt as seen in the two previous sections. Increasing the Lagrangian order 



© RAS, MNRAS 000, [T]-?? 



Non-linear evolution using LPT re-expansions 21 



and/or increasing the frequency of re-expansion reduces the differences. The right panel shows an anomaly, that the Cauchy 
differences for n = 2 and n — 3 lie on top of each other. This implies that convergence with Lagrangian order is broken. 

The lower panel compares momentum conservation error with (left) and without (right) frame shifts. The momentum, 
as defined in equation ( 64 1 , is calculated at initial and final times, P(to) and P/, respectively. The former is known exactly 
from the choice of initial conditions. After the evolution P/ is known numerically and it may be compared to the exact result 
P '/, exact = a(to)P(to)/df. We calculated the error as the absolute difference between the exact and approximate value summed 
in quadrature over the three directions. The figure (left) shows that conservation is achieved asymptotically with frame shifts. 
By contrast, there is no evidence for improvement when the frame shifts are omitted (right). Increasing the Lagrangian order 
and/or increasing the frequency of re-expansion reduces the error. 



3.7 Summary of convergence tests 

The tests indicate that the LPT re-expansion method works in principle and in practice for inhomogeneous initial conditions. 
This is the natural extension and generalization of our previous work on the top-hat model. These tests probe convergence 
when various combinations of grid size, Lagrangian order and/or re-expansion frequency are increased. Testing the rates 
separately has yielded results consistent with the underlying solution methodology (exponential in number of grid points for 
the spectral representation of spatial functions, exponential in Lagrangian order for the perturbation expansion, geometric 
in step size at fixed Lagrangian order). An essential part of the method is frame shifts. Ignoring them leads to violations of 
momentum conservation and breaks the theoretically expected convergence rates. 



4 CONCLUSION 

We have developed a new method to calculate the growth of large scale structure based on Lagrangian perturbation theory. It 
is designed to handle generic inhomogeneous perturbations on sub-horizon scales having non-relativistic velocities and weak 
gravitational fields. It can follow the system evolution until the formation of the first caustic. 

Many of the fundamental ideas were first discussed in the context of spherical top-hat cosmologies in NC11. There we 
demonstrated the existence of limitations on the convergence of the series expansion which is the basis of LPT. We showed 
that even if one could work to infinite order in the LPT expansion parameter there are bounds in time that limit how far 
forward a single step can extend. Such limitations are more restrictive than the onset of particle crossings. To circumvent the 
bounds we introduced the idea of re-expanding the solution in overlapping time domains ("LPT re-expansion"). Our work 
generalized LPT from a single-step to a multi-step scheme permitting evolution all the way up to the onset of multi-streamed 
flow. 

This paper extends the analysis from spherical top-hat cosmologies to generic initial conditions in the context of Newtonian 
cosmology. A new and important theoretical point is the existence of frame shifts. Frame shifts originate because the starting 
point of traditional derivations of the LPT series involves taking derivatives of the geodesic equation or its equivalent. The 
perturbative theory and solution obtained in this way is insensitive to time-dependent translations and will generally not 
satisfy the original equations of motion unless this freedom is fixed (no such freedom was present in the spherical cases studied 
in NC11). We introduce two conceptual frames to explore the situation. The OBS frame is the frame in which the physical 
initial conditions are specified and the answer sought. The EB frame is the frame in which the perturbative solution is obtained 
via the formal LPT expansion described in EB97. We describe in detail how to transform from one frame to the other. The 
relationship between OBS and EB frames is non-trivial. For example, it is not set by a physically conserved quantity. In many 
instances frame shifts are required, e.g. in a finite cosmology simulation, the EB frame will generically differ from that of the 
preferred OBS frame tied to the cosmic microwave background. At each step one must determine the relationship anew. We 
argued that ignoring the frame shifts generically leads to a breakdown in momentum conservation. 

The multi-step algorithm we have developed can be summarized as follows. Start with initial density and velocity defined 
in the OBS frame with Lagrangian coordinates equal to uniform Eulerian grid coordinates at time t. Transform quantities to 
the EB frame (the density remains the same, velocity differs by a constant). Use the EB formalism to advance the density and 
velocity in the EB frame from t to time t + St subject to the restriction imposed by the time of validity, the span over which 
the LPT expansion converges. Compute the frame shifts accumulated during the evolution and transform the EB solution 
back to the OBS frame at the end of the step. Interpolate the final OBS quantities onto the uniform Eulerian grid at time 
t + St. These quantities serve as the initial data with new initial Lagrangian coordinates equal to Eulerian coordinate for the 
next step. 

The algorithm was numerically implemented and tested. All tests were done in a fi m = 1 cosmology over time intervals 
that obeyed the predicted time of validity. Since these are meant to be 'proof of principle' tests we performed them on grids 
of small size. 

Some tests involve comparisons to exact answers. The solution for a simple, plane-parallel problem can be written down 
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without reference to LPT. Since first order LPT should produce an exact result this enables several different checks to be carried 
out. We showed that our first order LPT treatment with frame shifts matches the exact solution. We also established that 
ignoring the frame shifts gives rise to problems and inconsistencies. Several different lines of argument (involving momentum 
conservation, convergence, single versus multiple steps) substantiated the fundamental role of the frame shift even in the 
context of this simplest-of-all problems. Discrepancies introduced by ignoring frame shifts might go undetected for some 
specific initial conditions or in a single step scheme but they become clearly apparent even in this simple geometry for general 
initial data and/or while utilizing a multi-step scheme. 

As a second exact comparison we tested our results against the analytic example given in |Buchert et aT] ( p97| . This 
allows a check of the implementation of the higher order LPT terms. This provided a complete check of the non-zero (spatial 
and temporal) terms up to third order and a consistency check that the spatial parts of all the other terms vanish. The specific 
numbers of each type of term are given in the text. The frame shift, which vanishes for this particular example, is also checked 
by this test. 

Other tests involve situations where the exact answer is unknown and we focus on whether the convergence of the 
algorithm agreed with theoretical expectations. There are three control parameters: grid size N s , Lagrangian order n and the 
frequency of re-expansion N t . We first discuss theoretical convergence to the exact answer. Generally, N s and one (or both) 
of n and N t must increase simultaneously. The refinement must be made subject to the time of validity criterion. 

For diagnostic purposes we show it is possible to choose initial conditions and parameters such that errors are largely set 
by a single parameter. This allows us to study the nature of the convergence of that part of the calculation. To demonstrate 
convergence with grid size we represented the initial density and velocity with a spherical top-hat profile filtered with a 
Gaussian. We showed that for these smooth initial conditions, the errors due to the finite grid size decreased exponentially 
with N s , when other errors due to Lagrangian order were negligible. This was in accordance with the expected behaviour 
because spectral methods are used to solve the spatial equations. Likewise, when the grid size errors are very small the 
convergence with Lagrangian order is exponential in accordance with the basic understanding that LPT is a power series 
expansion in the perturbation parameter. We show that there is a sharp transition from convergence being controlled by grid 
size N 3 to it being controlled by Lagrangian order n. 

For the remaining studies, we set N a large enough that we could focus on the convergence with Lagrangian order n and/or 
the number of time steps Nt. We considered initial data based on a random Gaussian field. Two types of initial conditions 
were considered: Zel'dovich initial conditions, where the velocity was proportional to the acceleration field, and generic 
initial conditions, where the initial velocity includes both longitudinal and transverse components. Though unanticipated we 
discovered that the rate of convergence of these complicated examples scaled nearly identically as the rate of convergence for 
the simple top-hat cosmologies (whose rates were analytically derived and numerically confirmed in NCff). To summarize, 
the convergence with n is exponential and that with Nt is algebraic. The internal consistency for solutions achieved when 
varying n and Nt (a common asymptotic solution and Cauchy differences which diminish at the theoretically anticipated 
rate) provides strong evidence that the system is being solved correctly. We also demonstrated that when the frame shifts are 
omitted the convergence tests fail and conservation of momentum is violated. 

We explored how Kelvin's theorem is satisfied for situations in which the initial Eulerian vorticity vanishes. Kelvin's 
theorem implies that the Eulerian vorticity vanishes at all future times prior to the formation of caustics. Since this is an 
exact result it serves as a useful check on the end-to-end solution, i.e. the transformations between OBS and EB frames and 
the LPT calculation within the EB frame. We confirmed numerically that the magnitude of the Eulerian vorticity converges 
to zero in the asymptotic limit just as it should. This is a non-trivial check as can be seen from the following considerations. 
If the initial vorticity vanishes then vorticity as measured in the Lagrangian frame is generated at third and higher orders 
(for single step Lagrangian solution; BE93 and B94). What happens to these contributions in the OBS frame? They should 
vanish as the calculation order increases. The key is that vorticity as measured in the Lagrangian frame is not the same as 
vorticity as measured in the Eulerian frame. This can be seen in several ways. For example, a second order LPT calculation 
for the initial conditions described above yields zero vorticity in the Lagrangian frame (to second order) but this still gives 
rise to a non-zero vorticity in the Eulerian frame because of the transformation between the two frames. Alternatively, in a 
multi-step scheme, Eulerian vorticity is generated even with a first order LPT scheme. The appearance of Eulerian vorticity 
at a given finite order does not violate the Kelvin circulation theorem as long as the order of that contribution increases with 
the order of the calculation itself. This is the situation that is indicated by the numerical results. (Likewise, increasing the 
number of steps causes the magnitude of the Eulerian vorticity to diminish.) Prior studies have examined the density and/or 
the Lagrangian displacement and, to the best of our knowledge, these aspects of the solution have never been tested. The 
frame transformation is needed to achieve this consistency. 

We now consider the future utility of this new method. 

The choice of the exact number of time steps, Lagrangian order and grid size which govern the numerical errors will 
ultimately depend on the application at hand. The numerical errors can be made as small as desired. However, we have not 
addressed the important drawback of LPT, which is its inability to model physics beyond shell crossing. Real cosmological 
applications will be limited by the occurrence of multi-streaming. Any approximation to account for the velocity dispersion 
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( Adler fc Buchert|199"9}|Buchert, Dommguez, fc Perez-Mercader|1999||Morita fc Tatekawa|2001||Tatekawa|2005[) or alterations 



to the basic dynamics (such as the adhesion approximation; see reviews |Sahni fc Coles|1995| IGurbatov, Saichev, & Shandarin 
20121 will introduce physical errors as opposed to the numerical errors with which we have dealt. The algorithm presented 



in this paper will require additional development to deal directly with shell crossing for most cosmological applications. 
Alternatively, a good physical understanding of the nature of the errors introduced by various approximations may suffice. 
In general terms there are three main advantages of the LPT re-expansion method presented in this paper. (1) The smooth 



treatment of the density and velocity fields frees the results from particle-related shot-noise contamination (see Joyce, Marcos, 
& Baertschigcr 2009 for a detailed discussion of this issue in the context of numerical simulations). (2) The method can be 
fine-tuned to control the numerical error by changing the number of time steps and/or the Lagrangian order. (3) The scheme 
is designed to deal with general initial conditions including cases where the velocity fields may have a rotational component. 
Although this generality has always been implicit in LPT's formal development (Moutar de et al.||1991 B92), it has mostly 
been omitted in applications to real cosmological problems because vortical modes decay away in expanding cosmologies. 
Working in full generality should permit an LPT based scheme to handle arbitrary initial data which develops during the 
course of later phases of evolution. In addition, given the recent interest in using information from peculiar velocities to 
constrain cosmological parameters, it should prove useful to be able to handle these modes even in the quasi-linear regime. 

We envision several applications of our method, particularly where accurate results are desired on quasi-linear scales. We 
mention a few of them below. 

One important application is the Baryon Acoustic Oscillation (BAO) signal reconstruction. The detection of the BAO 



signal in the SDSS survey (Eisenstein et al. 2005[ ) implies that measurements of this sort will provide a powerful means 
to constrain cosmological parameters. Flows on quasi-linear scales cause the observed signal to have lower amplitude and 
peaks slightly shifted with respect to linear theory predictions. As observational probes improve, a sub-percent level precision 
determination of the BAO scale will require very accurate reconstruction ( Percival et al.|2010 l. Already techniques based on 



single step LPT have been proposed (e.g. Eisenstein et al. 2007 Matsubara 



2008 



Tassev & Zaldarriaga 20121. Our multi- 



step re-expansion algorithm can be used to achieve the needed level of error through control of a combination of time-step 
and Lagrangian order. Error control in perturbation techniques is of course an important element for the establishment of 
cosmological constraints ( |Carlson, White, fc Padmanabhan||2009[ ). 

Another important application is modelling the peculiar velocity field and investigating the density-velocity relation. The 
peculiar velocity information is encoded in the observed redshift space distortions (RSD) seen in galaxy surveys which recent 
measurements (e.g. Guzzo et al.|2008 l have demonstrated have the power to distinguish between various cosmological models. 
In linear theory, the density and velocity divergence are proportional and consequently, the ratio of the redshift space to real 
space power spectra is a constant. Numerical simulations show a breakdown of that prediction, starting at scales as large as 
k ~ 0.03/i -1 Mpc (Jennings, Baugh, & Pascoli 20111, which would normally be regarded as well within the linear regime. 



Understanding these scales is important. Analytical methods based on perturbation theory (for e.g., Bernardeau 1992; Kitaura 
et al.|[2012 l or based on spherical dynamics (Bilicki & Chodorowski 2008; Nadkarni-Ghosh 20121 give approximate answers 
but definitive answers will require a more simulation-based approach. LPT re-expansion will be ideal since it is intrinsically 
smoother than particle-based approaches. The current code can track all components of the velocity field enabling detailed 
investigations of RSD power spectrum as has recently been discussed by |Zhang, Pan, fc Zheng] ( |2012[ ). 

Other important applications include the evolution of non-Gaussianity and the estimation of the growth of voids. If 
and when shell crossing is accurately handled the method of LPT re-expansion may nicely complement existing methods for 
modeling the growth of large scale structure. 
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APPENDIX A: MATHEMATICAL TRANSFORMATIONS 



The section below outlines the transformations that are performed on equations (111 and (121 to change the dependent 
variable from r = yeb to X. For simplicity we drop the subscript 'EB'. The l.h.s. of equation (111 is 

V,.r=^ = |^. (A1) 
on dXi dri 

Einstein's repeated summation convention is followed throughout. The inverse transformation from X-space to r-space is 
given as 

dXi 1 



dr-j dr k 
~dr~ ~ 2j eimneijk dX m d~X~ n ' 



where 



J = Det 



dn 
dX; 



^abc 



dr\ dr2 drs 



drt dr p dr q 
dXjdXidXn 



and tijk is the usual Levi-Civita symbol. Define L[A, B, C] = eimqtij, 



OAj 9Bj dC k 
dX t dX m dX q 



. This gives 



_ .. _ 1 dr\ drj dr k _ 1 

V r -r - — e lmn e l]k — — — - — 

The divergence equation then becomes 



L[f,r,r] 



2J 



-4ttG 



p m ,oao 3 (l + S(X,t )) 
J 



+ Px,o(l + 3w) 



ao \3(i+») 



(?) 



Using the standard definitions of f2 m .o, flx.o, Hq, multiplying by J and noting that J = |L[r, r, r], the equation is 
L[r, r, r] = -3H 2 n m>0 a 3 (1 + 5(X, t )) - ^ (1 + 3w)Q x ,o (^) £ [ r > r, r]. 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



The gravitational field is irrotational in Eulerian space so (V r x r) = 0. Consider the component of equation ( 12 1 along 

dr k dr r dr 3 



the i-th direction. Again use ( A2 1 to write 
dr k 



(V r x r)i = e ijk 



dt 



^ijk 



Using the identity tij k 
dr k dri dr k 



dr k 8X t 
dXi drj 

dr k dr k dri 



2r ^ lmn e jr s dx Q Q - 0. 



SirSks — SisSkr in the above equation 



= 2ei r , 



dr k dr k dn 



= 0. 



(A7) 



(A8) 



" rm dXidx m dx n l " ul dx l dx m dx n dx l dx m dx n 

Here we have set each component of the vector V r x r in the r basis equal to zero. One can also express the components of 

dX " 

V r x f in the X basis where the two bases are related by fi = -g^-X p . This gives 

0. 



dr k dr k drj dX p - 



(A9) 



dXi dX m dX„ dn 

But dri/dX n ■ dX p /dri — S pn . Since the basis vectors are all independent the individual components must be zero. The 
simplified condition for each n then becomes 

dr k dr k 



= 0. 



dX t dX„ 

Defining T 9 [A,B] = ei mq g^ k it is rewritten as T[r, r] = 0. L and T are linear in each of their variables. 



(A10) 



APPENDIX B: COMPUTATIONAL DETAILS OF THE EB SCHEME 

In the EB frame, the equations to be solved have the following form 

A L [V x -p (1) ] = ~H o 2 {l m>o a 3 6(X,t o ) (Bl) 

A T [VxX P (1) ] = (B2) 

Dt [Vx ■ P (n) ] = S (n ' L) (B3) 

A T [V x xp<"»] = S<"' T \ (B4) 
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where n refers to orders higher than the first and 

'a N3(i +u ,) 



Dt = haa+^a 2 Ho 2 (l + 3w)Q x ,o(^) "+a 2 ^J (B5) 



Dl = + (B6) 

The explicit form for the source terms is 

1„ 3 TT 2/1 /a \ 3 ( 1 +-) 



d 2 



a,/3 

- ^ aL[pW,p(«X]- ^ i£[pW,pW jP W] 

a.-\-fi = n q: + /3+'7 — n 

- \H 2 (l + 3w)n x>0 (^y a+W) L[ p ("\p^,p^) (B7) 

ck,/3,7 

cn+/3-Hy— h 

S (n,T) = _ f[p (Q) ,p ( ' 9) ], (B8) 

a,/3 
a + — n 

where a, /}, 7 can take any values from 1 to n — 1 and they add up to n. These are subject to initial conditions 

at first order p (1 ' L/T) (X, to) = 0, p (1 - L/T) (X, to) = v i/T (X, to) (B9) 
at higher order p (n < L/T) (X,t ) = 0, p (n ' i/T) (X, t ) = 0. (BIO) 



Since the spatial and temporal operators in equations ( Bl I - (B4| commute, the displacement field can be written as a 
linear combination spatial vectors with purely time dependent coefficients 

Z n,L/T 

p (n,L/T) _ b {n,L/T) F (n,L/T) ^ 

i=l 

where b and F are functions of t and X respectively i.e., b = b(t) and F = F(X) and Z Ui l {Z ni T) denote the number of 
independent longitudinal (transverse) terms at the n-th order. The superscripts denote the order and type of the term. Using 
this decomposition, one can denote the higher order source terms as 

s (n,L) = J2 h\ n ' L) (t) • L[F a ,F^,F J ] (B12) 



i=l 



Y ^ n ' T) W-T,[F Q ,F' 9 ] (B13) 



respectively, where the superscripts of Fs take any value between to n - 1 (0 corresponding to F(X) = X) and add up 
to n. It is possible to calculate all the independent spatial source terms symbolically at all orders using the symmetries and 
properties of the L and T from which the numbers Z n ,L and Z n ,r can be determined and temporal source functions h(t) can 
be extracted. 

The initial perturbation is described by one scalar field 5(X, to) corresponding to the initial density field in the EB frame 
and two vector fields v L (X, to) and v T (X, to) corresponding respectively to the curl-free and divergence-less parts of the initial 
velocity field (in the EB frame). The initial acceleration vector can be constructed from the initial density field via Poisson's 
equation and is also curl-free. This gives Z\^l = 2 and Z\ t T = 1. Table [BT] shows the number of terms at various orders. If 
the initial conditions start with zero transverse velocity then Zi,t = 0. In this special case the total number of equations to 
solve reduces from 9 to 4 at second order and from 64 to 20 at third order. 

Substituting equation (Bill up to first order in equations ( Bl I and (B2| gives, 

D^' L) [V X -F^ L) ] + D^' L) [\7 X -F^ L) ] = -^H 2 Q mfi a 3 8(X,t ), (B14) 
D t T &i 1,T) [Vx x F^ T) ] = 0. (B15) 



The initial conditions that p' 1 ^ satisfies are given by equations (B9|. There is a choice to be made in how these initial 
conditions translate to conditions on the temporal and spatial functions. We choose to set them as shown below. 
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Table Bl. Number of transverse and longitudinal terms as a function of Lagrangian order and type of initial conditions. 



Order Zel'dovich initial conditions 


v T = 






v T ^0 


n — 1 Zjl — z, — U 


Zjh — 2, 


— U 






n=2 Z L = 3, Z T = 


= 3, Zt 


= 1 


y,. 


= 6, Z T = 3 


n=3 Z L = 10, Z T = 6 


Z L = 12, Z T 


= 8 


z L -- 


= 37, Z T = 27 


Vx-FS u) =i(X,i ) F^ L) 


= v L (X,* ) 


F 


(1,T) _ 
1 — 


v T (X,t ) 



&£ 1,Z0 (to) = 
b[ 1 ' L) (t ) = 



f-ffo fin 



,0«0 



Dfb. 





= 

= 1 



£> t T #' T) = 0. 
b[ UT) (t ) = 
^' T) (to) = l 



Substituting equation (Bll I in the higher order equations equations (B3l and (B4| and imposing equation (BIO I gives 

V7_. . v (n,L) = £.[ F O iF /3 iF 7] y x x F ^,T) = f i [F Q ,F' 9 ]. 



D^b { "' L \i) 
ti n ' L \to) = 



h. 



(u,L) 



(*) 



A T &f' T) (i) 



ft 



(n,T) 



(i). 



(*o) 







(to) = 0. 
(to) = 0. 



All the spatial solutions are subject to periodic boundary conditions and are solved using Fourier transforms. 



APPENDIX C: FORMAL REQUIREMENTS OF THE EB SOLUTION 

The formal derivation by Ehlers and Buchert (EB97) demands that the displacement in the EB frame satisfies 
p(X,t)d 3 X = 0. 



(CI) 



for all times t. At the initial time the definition of the Lagrangian coordinate implies p(X, to) = and equation (CI I is 
trivially satisfied. We want to prove that if the initial data satisfies 



{8(X,t ))x = 
(v(X,to))x = 



(C2) 
(C3) 



and the system is periodic, then equation (CI I holds at all times. We refer to these as the 'zero-mean' initial conditions. 



Proof. The displacement p(X,t) is a time-dependent linear combination of some spatial functions F(X). We will prove 
that each F(X) satisfies (F(X))x = and therefore (p(X,t))x = 0. 

At first order, the spatial function is either a solution of Vx • F(X) = <5(X, to) or F(X) = v(X,to)- Since the initial 
density integrates to zero over the volume and the system is periodic, the solution also has the same property. Therefore, at 
first order all spatial Fs has mean zero. 

At higher orders the Fs are periodic solutions of Vx ■ F = 5 ,L (X) or Vx x F = S T (X), where S , S T are combinations 
of lower order Fs. So if S L (X) and S T (X) have zero mean then periodicity will imply that each higher order F has mean 
zero. Here we present the proof only for the longitudinal source terms S L (X.) = L[F a , F' 3 , F 7 ]. The proof for the transverse 
sources T[F a ,F' 3 ] follows a similar calculation. 

Let I — (L[F a , F* 3 , F 7 ]}. To keep in mind the periodic nature of the system, we refer to space as the 3-torus T 3 . Write 
I = yttjkl'ijk, where 



rl 



cfX 



d F? BFf dFl 
dX p dX q 8X r 



(C4) 



Here the spatial functions F can either be periodic functions of X or equal to X. Since these are higher order source terms 
and at least one of Fs is not X. Without loss of generality let this correspond to the term dFk/dX r . Integrating by parts over 
X r and using periodicity gives 



ji 




cfX 



d 2 F,° 



dx r dx„ / ox, 



dF: 



d 2 F? 
dX r dX q 




(C5) 
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The first integrand is symmetric under the exchange r «-> p and the second under the exchange r «-> q. e P q r is antisymmetric 
under these exchanges. Therefore, = and hence I — 0. Simiiarly one can show that the source terms for solving the 
transverse vectors also integrate to zero. Thus equation (CI I holds for all t. 



APPENDIX D: SETTING THE INITIAL CONDITIONS ALONG THE ZEL'DOVICH CURVE 



It is generally common in cosmology to assume that the initial velocity is proportional to the initial acceleration ( Zel'dovich 



|1970[ ); the proportionality constant is set by requiring that there be no perturbations at the big bang. This is also equivalent 
to having no decaying modes i.e. no negative powers of t in the displacement. At very early times (recombination) the universe 
is purely matter dominated with Q = 1 and the first order displacement evolves as (B92) 



J 1 ) 



(I) 



-1/3 



2a _(i,£) 



3io l 



t \ 4 / 3 r 



3ao F (i,£) 



3to L 

If" 



2/3 



■ao — 



F (l,i)_ 



Setting the coefficient of the negative power of t to zero gives 



a 0„(l,L) 



3 t 



F^ J = -a(t )F\ 



(Dl) 



(D2) 



Note that this does not guarantee that there are no decaying modes from the higher order solution. The terms (i/io) -1 ^ 3 and 
(t/to) 4 ^ 3 arise from the homogeneous part of the solution and each higher order solution will contain them since the temporal 
derivative operator is the same at all orders. An alternate way is to choose the initial velocity such that the scaled velocity 
divergence 5v = (3Ho)~ 1 V roBS ■ vobs and 5{X.,to) at each point satisfy the non-linear Zel'dovich relation based on the 
top-hat analysis (NC11). In this case the relationship in equation (D2l is not satisfied; decaying modes will be present even 



at first order. The former prescription is used more often in literature and in this paper. We also checked that the difference 
in the two ways of setting the initial conditions is very small (~ 10 -8 ) when the starting epoch is z ~ 1000 (recombination). 



APPENDIX E: ID PERTURBATIONS 

This appendix proves the results for the 1-D case stated in §3.3| 

El Frame shift for a single step 

Since the problem is essentially 1-d, we will denote all spatial functions as scalar functions. The Lagrangian coordinate in the 
observer (EB) frame is denoted as Y {X); Y = X. The initial density and velocity in the observer frame are 5oBs(Y,to) and 
vobs (Y, to) respectively. The initial data in the EB frame is 

Seb(X, t ) = 5 Bs{Y,to); v EB (X, t ) = v bs(Y, t ) - « c ,o, (El) 

where u c ,o = {vobs(Y, to)). The acceleration in the observer's frame is given by Fg(Y,to) where dFs/dY = 5obs(Y, to). The 
densities are the same, hence F$(X,to) = Fg(Y,to). Hereafter we drop the subscripts on the density since it is the same in 
both frames. 

The displacement in the EB frame at any later time is given as 
p(X,t) = b 1 (t)F 1 (X) + b 2 (t)F 2 (X), (E2) 
where F\(X) = Fg(X,to) and F 2 (X) — VEB{X,to). The physical coordinate 

r EB {X,t) =aX+p(X,t) (E3) 
and the corresponding Jacobian is 



The comoving coordinate is xeb = teb/c. Restricting to ID and substituting equation (E3l and (E4i in equation (551, the 
frameshift equation is 



d ( o.dAx\ 1 
dt 

Since p is a periodic function of X , this gives 

d ( tdk.x\ 1 
dt 
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Substituting for p from equation ( E2 1 

-C 

dt \ 



,dAx 



dt J 



L 



(b 1 b 2 F 1 Fi + b 2 b 1 F 2 Fl)dX, 



(E7) 



where prime denotes differentiation w.r.t. X or Y . Using the equations for 61 and 62 given in Appendix |B] applying the 
product rule for derivatives and using periodicity of the spatial functions gives 



( = ~ ( 3a+ 2 (1 + 3u,)anx 'Hi) \b 2 - (F^^dX (x F s (Y,t )v' OBS (Y,t )dY, (E8) 

with initial conditions Ax(to) = and Ax(to) = v c ,o/ao- Consider the case where voBs(Y,to) oc Fs(Y,to). The condition that 



the mean density be zero implies {Fg(Y, to)) = « c ,o = 0. Thus the initial conditions as well as source term of equation (E8l 
are zero and hence Ax = for all t. Thus, single step LPT with Zel'dovich initial conditions requires no shifts. 



E2 Maintaining the Zel'dovich condition 

Let the initial time at the start at the first step be to and consider taking a second step at time t\. Let Yo, Y\ be the Lagrangian 
coordinates in the observer frame at time to, ti respectively. Y± is related to Yq as 

Yi=Y + a - 1 p(X,ti) + Ax(ti) (E9) 
If VoBs(Yo,to) oc Fs(Yo,to), then v c ,o = and VBB(X,to) = VoBs(Yo,to)- Combining with equation ( |E2[ ), this gives 
p(X,ti) oc Fs(Y ,to) and p(X,ti) oc F s (Y ,t ). (E10) 
The velocity in the observer frame at the start of the second step is 



VOBs(Yl,ti) =p(X,ti) 



d(ti) 
a(ti) 



o(X,ti) + aAx(ti). 



The acceleration at the start of the second step is 
Fg(Yi,ti) = / o(Yi,ti)dYi = jJxT) 



where J(X, t) is given by equation |E4|. Transforming from Y\ to X(= Yq) using equation (E9l 
Fi(¥i,ti) = J(l + d(Y ,to))dY - I (} + \%) dX 
= Fs (Y ,to)- P -^ 1 . 



(Ell) 



(E12) 



(E13) 



(E14) 



When the initial acceleration and velocity are proportional, the above equation can be combined with equation (ElOl to give 



Fg(Yi,ti) oc Fs(Yo, to). Since the frame shifts are zero, equation (ElOl and ( Ell I give voBs{Yi,ti) oc Fs(Yo,to). This proves 
that if the initial acceleration and velocity are proportional they remain proportional at all subsequent times. 



E3 Conservation of mean velocity 

Using the periodicity of p(X,t), equation (E6l can be re- written as 

d ( 2 dAx\ d ( 1 f .dp 
di{ a ^t-) = dt [-L J P dX dX 

This gives 

; M «c,oao 1 (l f .dp 
Ax (t) = —2— -^[l I P^ dX 



dX 



The velocity in the observer frame at any time is given by equation I Ell I. The mean at any time t is 



{voBs(Yi,t)) = L I r ( , D ,( V,. t)d) 1 



p--p + aAx) 



adXl 



aAx+1 a\L h>JX- dX 



(E15) 

(E16) 

(E17) 
(E18) 
(E19) 



Substituting for Aa; from equation (E16l 
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(v Bs(Y,,t)) = V -^. (E20) 

Thus, for a 1-d system the mean velocity, like the mean momentum, is conserved modulo a '1/a decay' due to the Hubble 
drag. 



APPENDIX F: SPHERICAL TOP-HAT 

This section describes the details involved in setting up the compensated top-hat configuration for section [375] 

The exact compensated top- hat function consists of a spherical overdense region surrounded by a compensating underdense 
vacuum region. Let a and b be the radii of the overdense and compensating regions respectively. The initial density profile is 
given by 

{So sC |X| < a 
-1 asC|X|<6 (Ff) 
|X|>6 

We chose So = 10, a = 1/4 and b = ll 1//3 /4. The box length Lt ox is chosen to be two units in length centred around the point 
(1/10, —1/11, l/(27r)). The choice of parameters ensures that the entire profile is well represented within the box and the 
offset ensures that no special symmetry is exploited in the test. The profile is discontinuous at X = a and X — b. The Fourier 
transform of a discontinuous function has power at all wave numbers and the Gibbs phenomenon prevents such functions 
from being completely represented by Fourier transforms on any finite size grid. The initial profile is hence smoothed with a 
Gaussian filter of width a to give 

= jT" s(x,t o) ■ %Srf ) d3x - ^ 

The smoothing is performed analytically in real space and the resulting function is numerically evaluated on the grid. To 
ensure periodicity of the initial data, contribution of the 26 nearest neighbours cells is added. The contribution of cells beyond 



the nearest neighbours was zero to machine precision. For the density profile in S3.5 the smoothing parameter is a = 1/12 



APPENDIX G: GENERATION OF THE TRANSVERSE VELOCITY 

The Eulerian vorticity in the observer frame is V roBS x vobs- 

_ dvoBS,k ,„.s 

Vr OBS X VOBS = tijk-F, • (,<jf) 

OroBSJ 

LPT expresses vobs as functions of the Lagrangian variable Y and not vobs- Using the relations in Appendix [A] and 
properties of the Levi-Civita tensor gives 

(^r OBS X VOBs)i = ^-j£ijk£jpq£lmnrp,mr q ,nVk,l (G2) 

= Tj (SkA, ~ M* 9 )ejmnr P , m r g ,nt;w 

= -^-j£imn(rk, m ri, n Vk,i — r itm rk,nVk,i) (G3) 

= -j(-m,ln{ri^ m rk,nVk,l), (G4) 

where the 'comma' denotes differentiation with respect to the Lagrangian variable Y = X. Substituting for vobs and 



voss from equations (591 and (611 gives rise to three types of terms in the expression: {eukPk,i, tukPk.i}, {tiinPk,nPk,i} and 
{timnPi,mPk,nPk,i\ ■ The first kind corresponds to the vorticity in the Lagrangian frame while the remaining two have no 
obvious physical interpretation. In general, all three types of terms are non-zero. For a first order LPT calculation (i.e., the 
displacement is accurate to first order in the expansion parameter and leading errors are second order), p is always proportional 
to p and the anti-symmetry of the Levi-Civita tensor ensures that all three types of terms vanish. The second order LPT series 
does not have this feature. In particular, the second order displacement term in the series p' 2 ' is not proportional to the first 
order p' 1 ^ and hence combinations of these (which are third order) will in general not be zero. Thus, the Eulerian vorticity 
will in general have non-zero contributions which are third order and higher (up to order 6) in the expansion parameter. Thus, 
when calculations are performed with second order LPT, the Eulerian vorticity is non-zero although the Lagrangian vorticity 



© RAS, MNRAS 000, [T]-?? 



